Since a mixed state is a normal statistical combination of many pure states, is it possible to find the list of pure state that makes up the mixed state?
I realise there are many different pure states that can combine into the same mixed state, so "disassembling" may produce infinitely many solutions. Is it possible to perhaps find one of the solutions? Or the "most probable" pure state under the classical understanding?
I realise there are many different pure states that can combine into the same mixed state, so "disassembling" may produce infinitely many solutions.
This is indeed true. The decomposition of a density matrix into a classical distribution is highly ambiguous (as I explained at depth in this previous thread).
However, the loose idea of the "most probable" pure state that helps constitute a given density matrix can indeed be given a unique, rigorous understanding. In essence, a density matrix $\rho$ is a self-adjoint, positive-semidefinite, trace-class operator such that $\mathrm{Tr}(\rho)=1$. The first of those qualifiers is the key, because as a self-adjoint operator you are guaranteed* an eigenvector decomposition of the form $$ \rho = \sum_n p_n |\phi_n \rangle\langle \phi_n|.$$ There the $|\phi_n\rangle$ are guaranteed to be orthonormal, and they are unique (up to degeneracies). Moreover, the probabilities $p_n$ satisfy $p_n\geq 0$ and $\sum_n p_n =1$.
This structure then lets you take the $|\phi_n\rangle$ with maximal $p_n$ as the "most probable" pure state within the probability distribution.
Moreover, because that maximal $p_n$ coincides with the operator norm of $\rho$, that "most probable" state can also be defined as the state $|\phi\rangle$ such that $$ \langle\phi|\rho|\phi\rangle $$ is maximal over the entire Hilbert space (as can be shown by expanding this $|\phi\rangle$ over the eigenvector basis and then calculating this expectation value).
*modulo mathematical complications in infinite dimensions which require more work, but the assumptions are very strong (the operator is bound, and likely compact) so you get a strong spectral theorem in return.
One possible decomposition is given by the eigenvalue decomposition $$ \rho = \sum p_k \vert\psi_k\rangle\langle\psi_k\vert\ . $$ Any other decomposition $$ \rho = \sum q_\ell \vert\phi_\ell\rangle\langle\phi_\ell\vert$$ is related to this one through $$ \sqrt{q_\ell}\vert\phi_\ell\rangle = \sum v_{\ell k}\sqrt{p_k}\vert\psi_k\rangle\ , $$ where $(v_{\ell k})$ is an isometry, i.e. $$ \sum_\ell \bar v_{\ell k'}v_{\ell k} = \delta_{k,k'}\ . $$ (Note that this is an if and only if, that is, a full characterization of all ensemble decompositions of $\rho$.)
