How to compute the pure extensions of a given mixed state? Let us consider any pure state $|\psi\rangle\in\mathbb{C}^n\otimes \mathbb{C}^n\otimes \mathbb{C}^n$. Its reduced bipartite density matrix represent a pure state or mixed state depending on whether $|\psi\rangle$ is entangled or not (exactly how it is entangled, on which system we take the partial trace, etc). 
My question is given any arbitrary (mixed) state $\rho\in\mathcal{B}(\mathbb{C}^n\otimes \mathbb{C}^n)$, can we find out a  pure state $|\psi\rangle\in\mathbb{C}^n\otimes \mathbb{C}^n\otimes \mathbb{C}^n$ (or in some suitable higher dimension which needs to be determined) such that $\rho$ is the reduced density matrix of $|\psi\rangle$. In particular, I do not want only an existential result,  I also want to an algorithmic method to determine such  $|\psi\rangle$. Obviously such state will not be unique. Advanced thanks for any help. 
 A: *

*The fact that every mixed state $\rho$ acting on a finite dimensional Hilbert spaces can be viewed as the reduced state of some pure state $|\psi\rangle$ on a bigger Hilbert space is known as purification, see this Wikipedia page, where also the algorithm is given. 

*In OP's case of $$\rho~\in~\mathcal{B}(\mathbb{C}^n\otimes \mathbb{C}^n),$$ one may choose a pure state $|\psi\rangle$ in the following Hilbert space
$$|\psi\rangle~\in~\mathbb{C}^n\otimes \mathbb{C}^n\otimes \mathbb{C}^n\otimes \mathbb{C}^n.$$
A: It is a general result that any mixed quantum state can be viewed as the reduced state of a pure state in a larger dimensional Hilbert space. This is referred to as a purification, and some people even refer to the power of this idea as the "Church of the Larger Hilbert Space". 
There exists a canonical way of constructing a purification, which has the advantage that it immediately demonstrates what the minimum dimension of the larger Hilbert space must be. So here is how you construct a purification: Let $\rho$ be a state in a Hilbert space of dimension $d$ with eigenvectors $\{|\phi_i\rangle\}_{i=1}^d$. Then 
$|\psi\rangle=\sum_{i=1}^d|i\rangle\otimes|\phi_i\rangle$
is a purification of $\rho$, where $\{|i\rangle\}$ is an orthonormal set of vectors. 
Therefore, you can see that it is always possible to construct a purification and moreover, the enlarged Hilbert space in general must have dimension at least the twice as big as the original one.
A: You correctly say that such a state won't be unique. Indeed, for example, pick the 50%-50% statistical mixture of two random pure basis states in ${\mathbb C}^n\otimes {\mathbb C}^n$. Each of these two states may be entangled with one of two orthogonal states in the third ${\mathbb C}^n$; however, what these two states in the third factor Hilbert space happen to be is completely undetermined.
So even if you had a constructive method to find the pure state in the three-part Hilbert space, it would fail to generate unique results.
However, there is another, much more serious problem with your proposal: in almost all cases, it has no solutions at all. Indeed, it's easy to demonstrate this fact by a simple counting of degrees of freedom. Pure states in ${\mathbb C}^n\otimes {\mathbb C}^n$ are specified by $n^2$ different complex numbers (one of them, the overall complex normalization, is unphysical).
Similarly, a general density matrix on this space is a $n^2\times n^2$ Hermitian matrix so it contains $n^4$ independent real parameters (one of them is the trace which should probably be set to one). However, that's larger than $2n^3$ number of real parameters coming from $n^3$ complex parameters of a wave function in ${\mathbb C}^n\otimes {\mathbb C}^n\otimes{\mathbb C}^n.$ At least for $n\gt 2$, it's larger. So up to a measure-zero subset of cases, you won't be able to find any pure state that reduces to the given mixed state. The diversity of the required results (density matrices) is much larger than the diversity of the ingredients (pure three-block states) that you may use to produce the desired outcome.
Of course, if you only had a density matrix for one of the three blocks, and not two of them, you would be able to solve it. At least the counting of the parameters wouldn't make the existence of a solution for a generic density matrix impossible.
A: Given an arbitrary state $\rho\in\mathcal B(\mathbb C^n)$ whose eigendecomposition reads
$$\rho=\sum_{k=1}^{\mathrm{rank}(\rho)} p_k|\psi_k\rangle\!\langle \psi_k|,\quad p_k>0,\tag A$$
The set of its purifications is the set of vector states $|\Psi\rangle\in\mathbb C^{n}\otimes\mathbb C^{m}$ for $m\ge \mathrm{rank}(\rho)$ whose SVD (called Schmidt in this context) decomposition has the form
$$|\Psi\rangle=\sum_{k=1}^{\mathrm{rank}(\rho)}\sqrt{p_k}\,|\psi_k\rangle\otimes|u_k\rangle,\tag B$$
where $\{|u_k\rangle\}_k\subset\mathbb C^m$ is an arbitrary orthonormal set of $\mathrm{rank}(\rho)$ elements of $\mathbb C^m$.
From this we can conclude a number of things:


*

*A generic state $\rho$ can have rank up to $n$, and so the minimal dimension of the purification space needed to accommodate the purifications of arbitrary states is $m=n$.

*The bipartite structure in the state used in the OP is not important in this discussion. If $\rho\in\mathcal B(\mathbb C^n\otimes\mathbb C^n)$, that just means that the dimension of the space is $n^2$, and everything else follows as shown in the general case in which we don't reference the partite structure of the states.

*Eq. (B) makes it very clear what the possible purifications of $\rho$ are: the freedom is all and only in the choice of an orthonormal set of $\mathrm{rank}(\rho)$ vectors from some arbitrary ancillary space (with the only caveat that this space needs to be large enough to accommodate this number of orthogonal elements).


As a final remark, note that what in this context is referred to as purification, from a mathematical point of view is equivalent to the characterisation of positive operators $B$ as those operators such that $B=A A^\dagger$ for some $A$.
That is, the problem of purifying a given state $\rho$ is akin to that of finding $A$ such that $A^\dagger A=B$ for some given $B\ge0$.
To see this, we just need to realise that the partial trace operation on a rank-1 projector, $\operatorname{Tr}_2(v v^*)$, is equivalently stated as matrix multiplication of the operators having $v,v^*$ as vectorisation.
More precisely, by this I mean that, for any (possibly rectangular) matrix $A$,
$$A A^\dagger = \operatorname{Tr}_2[\operatorname{vec}(A)\operatorname{vec}(A)^*],$$
where $\mathrm{vec}(A)_{ij}\equiv A_{ij}$ and $\mathrm{vec}(A)\in\mathbb C^{s\cdot r}$ if $A$ is an $s\times r$ matrix.
The connection with the bra-ket formalism and this is the identification $|\Psi\rangle=\operatorname{vec}(A)$.
