What is the minimum number of separable pure states needed to decompose arbitrary separable states? Consider a separable state $\rho$ living in a tensor product space $\mathcal H\otimes\mathcal H'$, with $\mathcal H$ and $\mathcal H'$ of dimensions $D$ and $D'$, respectively.
If $\rho$ is separable, then it is by definition possible to write it as a convex combination of (projectors over) separable pure states.
$\newcommand{\ketbra}[1]{\lvert #1\rangle\!\langle #1\rvert}$Because a state is Hermitian and positive by definition, we can trivially always write it in terms of its eigenvectors and eigenvalues as
$$ \rho = \sum_{k=1}^{D^{} D^\prime} \lambda_k \ketbra{\psi_k}, \quad p_k\ge0, $$
where $\rho|\psi_k\rangle=\lambda_k|\psi_k\rangle$.
However, $|\psi_k\rangle$ will in general be non-separable states.
What I am looking for is the decomposition of $\rho$ in terms of only separable states.
For example, a trivial case is $\rho=I/DD'$, which is easily seen to be decomposable as
$$\frac{1}{DD'}I=\frac{1}{DD'}\sum_{k=1}^D\sum_{\ell=1}^{D'}\ketbra{k,\ell}.$$
This shows that, to decompose an unknown state $\rho$ in terms of separable states, at least $DD'$ elements are required.
Is this number sufficient for any separable $\rho$?
In other words, what I'm looking for is the smallest $M$ such that a representation of the form
$$\rho = \sum_{j=1}^M p_j \,\ketbra{\alpha_j}\otimes\ketbra{\beta_j}$$
holds for all separable $\rho$.
More formally, this amounts to finding
$$\min\left\{M\in\mathbb N\,:\,\,\forall\rho\exists\{p_k\}_k,\{|\alpha_k\rangle\}_k,\{|\beta_k\rangle\}\,:\,\rho=\sum_{j=1}^M p_j \,\ketbra{\alpha_j}\otimes\ketbra{\beta_j}\right\}.$$
 A: First of all, your problem is a special version of a more general problem, namely finding the minimum number of states which minimize the entanglement of formation, this is, given a state $\rho$ on AB$\equiv \mathbb C^D\otimes \mathbb C^{D'}$, find the decomposition
$$
\rho = \sum_{i=1}^m p_i |\psi_i\rangle\langle\psi_i|
$$
which minimizes $\sum_i p_i E(|\psi_i\rangle)$, where $E(|\psi_i\rangle) = S(\mathrm{tr}_B(|\psi_i\rangle\langle\psi_i|)$, and find the minimum $m$ for which such a decomposition exists. 
Your problem is just the variant of this where the state has entanglement of formation zero.
This is a well-studied problem and in turn a special case of a so-called "convex roof construction".  Uhlmann, for instance, states that for any such problem, at most $(DD')^2+1$ states are needed for the optimal decomposition (Proposition 2.1).

It is likely that better bounds exist for the special problem of entanglement of formation, or the given problem of separable states.  I was unable to find any in the literature, but one should be able to prove one along the following lines:


*

*First, note that one can relax the optimization to all decompositions 
$$\rho=\sum p_i\rho_i\,\tag{1}$$ where one minimizes $\sum p_i S(\mathrm{tr}_B\rho_i)$, since the entropy is concave, i.e. the minimium will always be (also) attained on pure $\rho_i$. 

*Thus, we can instead consider decompositions of the reduced density matrix $\rho^A = \sum p_i \rho_i^A$ -- any such decomposition arises from a decomposition (1) of $\rho$ (e.g. by writing $p_i\rho_i^A$ as $M_k\rho M_k^\dagger$ with a POVM $M_k$ and applying $M_k\otimes I$ to $\rho$).

*Now consider an optimal decomposition $\rho^A = \sum p_i \rho_i^A$. If it has more than $D^2$ terms, the $\rho_i^A$  must be linearly dependent.  Thus, we can decrease the weight of some $\rho_j^A$ all the way down to zero by shifting the weights of all the other $\rho_i^A$ (keeping $p_i\ge0$!). Again, due to concavity, this will not change the average entanglement.

*We are now left with an optimal decomposition $\rho^A=\sum p_i\rho^A_i$ with $D^2$ terms.  This yields a decomposition of $\rho$, $\rho=\sum p_i \rho_i$, which minimizes $\sum p_i S(\rho_i^A)$ (as described in 2.). We can now decompose each $\rho_i$ in their eigenbasis (which has at most $DD'$ terms), which yields a total of $D^3D'$ terms. 

*There is likely space for improvement: For instance, one could rewrite each of the $\rho_i^A$ in a basis of pure states $|\phi_{k,i}\rangle\langle\phi_{k,i}|$. Such a basis has size at most $D^2+1$ ($D^2$ being the dimension of the convex space), and the coefficients are $\mathrm{tr}(\rho_i^A|\phi_k\rangle\langle\phi_k|)$ and thus positive. Again, convexity yields an optimal decomposition with pure $\rho_i^A$ and $D^2$ terms. It only remains to decompose the corresponding $\rho_i^B$, which results in a total of $(D^2+1)D'$ terms.
A: As stated above, Caratheodory's theorem provides an upper bound of $ (D D')^2 $, but it is not always strict. For example, for two qubits one never requires more than $ M=4 $ terms (see here). It turns out the general case is nontrivial, and has been studied in the literature. See for example this paper and the works cited there.
Personally, I am interested in the similar problem of determining the minimal number of necessary product terms, but where the single-party states are not required to be pure. This problem is less studied and also seems to be nontrivial (I had convinced myself in the past that $ \min \{ D^2, D'^2 \} $ terms should always suffice, but apparently I relied on a false argument).
A: It's obviously not always possible -- consider, for example, a pure entangled state. In this special case, the decomposition is unique and contains only one term, the entangled state itself $\rho=\lvert\psi\rangle\langle\psi\rvert$.
Given an arbitrary state, it may be tricky to determine if you can find a decomposition involving only separable states, but here's one suggestion that might help: you could check the eigenvalues of $\rho$.


*

*Any decomposition $\sum_j \lvert\psi_j\rangle\langle\psi_j\rvert$ can contain only $\lvert\psi_j\rangle$ orthogonal to all the eigenstates with eigenvalues equal to zero.

*In case any subset of the non-zero eigenvalues are degenerate and correspond to entangled states, you can try to construct linear combinations of them that are not entangled.
Here's an example: consider the state $\rho=\frac{1}{2}(\lvert00\rangle+\lvert11\rangle)(\langle.\rvert)+\frac{1}{2}(\lvert00\rangle-\lvert11\rangle)(\langle.\rvert)$. It has two distinct eigenvalues, $\frac{1}{2}$ and $0$, both degenerate. Therefore, no decomposition will include the terms $(\lvert01\rangle \pm \lvert10\rangle)(\langle . \rvert)$. Furthermore, since the eigenvalues for $(\lvert00\rangle \pm \lvert11\rangle)(\langle . \rvert)$ are degenerate, the density matrix will also be diagonal in any basis that uses a linear combination of these states. For example, we could use
$\{(\lvert 00\rangle+\lvert 11\rangle)\pm(\lvert 00\rangle+\lvert 11\rangle)\}= \{\lvert00\rangle,\lvert11\rangle\}$,
so the state can also be written $\rho=\frac{1}{2}\lvert00\rangle\langle00\rvert+\frac{1}{2}\lvert11\rangle\langle11\rvert$.
