# How to find all the possible density matrices $\rho_{AB}$ that produce a given pair ($\rho_A$, $\rho_B$) of reduced density matrices?

For a bipartite density matrix ($\rho_{AB}$), each of reduced density matrices ($\rho_A$ and $\rho_B$) are uniquely defined. But the converse is not true. That is, more than one $\rho_{AB}$ can have same $\rho_A$ and $\rho_B$. e.g- all Bell states have same reduced density matrix. Given $\rho_A$ and $\rho_B$, is there any systematic way to find all possible $\rho_{AB}$? Do all these states have any connection?

Given $\rho_A$ and $\rho_B$, choose bases $\{|n,k\rangle_A\}$ for $\mathcal H_A$ and $\{|m,l\rangle_B\}$ for $\mathcal H_B$, so that $$\rho_A = \sum_{n,k} p_n |n,k\rangle\langle n,k|_A,$$ i.e. $k$ iterates over the possibly degenerate subspaces. Similarly, for $B$, write $$\rho_B = \sum_{m,l} q_m |m,l\rangle\langle m,l|_B.$$ Since these are complete bases, then we can write $$\rho_{AB} = \sum_{n,k}\sum_{n',k'} \sum_{m,l}\sum_{m',l'} r_{nkn'k'mlm'l'} |n,k\rangle\langle n',k'|_A \otimes |m,l\rangle\langle m',l'|_B.$$

Finally, take the partial trace of this density matrix to get $$\mathrm{Tr}_B(\rho_{AB}) =\sum_{n,k}\sum_{n',k'} \sum_{m,l}r_{nkn'k'mlml} |n,k\rangle\langle n',k'|_A =\sum_{n,k} p_n |n,k\rangle\langle n,k|_A =\rho_A,$$ and similarly for $B$. That then gives you the requirement that the off-block-diagonal elements vanish, i.e. $r_{nkn'k'mlm'l'}=0$ under either of $n\neq n'$ and $m\neq m'$, but you can probably still have nonzero off-diagonal elements within each block. In addition, you also require the sums to check out: \begin{align} \sum_{k,k'} \sum_{m,l} r_{nknk'mlml} = p_n\\ \sum_{n,k} \sum_{l,l'} r_{nknkmlml'} = q_m. \end{align}

There's probably very little else that you can say, in general. Working in the natural-orbital basis (so you're matching against diagonal instead of general reduced tensity matrices) can eliminate one sum out of the game, but that's likely to be as far as you can go.

The brute force way of doing it:

1. Find a basis for your density matrix space. For a two-level/spin system the Pauli matrices may be convenient.
2. Expand the density matrix in this basis.
3. Compute the required reduced density matrices $\rho_A$ and $\rho_B$.
4. Impose $\rho_A = \mathrm{\{desired\quad reduced \quad matrix\}}$ and $\rho_B = \mathrm{\{desired\quad reduced \quad matrix\}}$ as restrictions on the basis expansion coefficients.
5. Solve this linear system of equations for your choice of parameters to be eliminated.
6. The result is a density matrix in terms of a number 1 of independent parameters.

1 In general, the number depends on which dimension your density matrix is and what your desired reduced density matrices are.