Considering a state that describes two subsystems $X$ and $Y$:
$$ \left|\psi \right\rangle =\sum_{i,k} \beta_{ik} \left| \theta_i^X \right\rangle \otimes \left| \lambda_i^Y \right\rangle$$
I am trying to show that it is possible to choose a basis for the space of states of $X$ and $Y$ such that the density matrix would reduce to the following form:
$$ \left| \psi \right\rangle \left\langle\psi \right| = \sum_{i,j} \alpha_{ij} \left| \phi_i^X \right\rangle \otimes \left| \psi_i^Y \right\rangle \left\langle\phi_j^X \right| \otimes \left\langle\psi_j^Y \right| \hspace{0,7cm}(1)$$
My attempt consisted in plugging in four identity operators that come from two new orthonormal basis: $\sum_i\left| \phi_i^X \right\rangle \left\langle\phi_i^X \right|= \hat{\mathbb{I}}$ and $\sum_j\left| \psi_j^Y \right\rangle \left\langle\psi_j^Y \right|= \hat{\mathbb{I}}$ such that:
$$ \left| \psi \right\rangle \left\langle\psi \right| = \sum_{i,j} \sum_{k,l} \beta_{ij} \beta_{kl}^* \hat{\mathbb{I}}\left| \theta_i^X \right\rangle \otimes \hat{\mathbb{I}}\left| \lambda_j^Y \right\rangle \left\langle\theta_k^X \right|\hat{\mathbb{I}} \otimes \left\langle\lambda_l^Y \right|\hat{\mathbb{I}}$$
$$ \left| \psi \right\rangle \left\langle\psi \right| = \sum_{i,j} \sum_{k,l} \beta_{ij} \beta_{kl}^* \left( \sum_n\left| \phi_n^X \right\rangle \left\langle\phi_n^X \right| \right)\left| \theta_i^X \right\rangle \otimes \left( \sum_{\alpha}\left| \psi_{\alpha}^Y \right\rangle \left\langle\psi_{\alpha}^Y \right| \right)\left| \lambda_j^Y \right\rangle \left\langle\theta_k^X \right| \left( \sum_n\left| \phi_n^X \right\rangle \left\langle\phi_n^X \right| \right) \otimes \left\langle\lambda_l^Y \right|\left( \sum_{\alpha}\left| \psi_{\alpha}^Y \right\rangle \left\langle\psi_{\alpha}^Y \right| \right)$$
$$ \left| \psi \right\rangle \left\langle\psi \right| = \sum_{i,j} \sum_{k,l} \beta_{ij} \beta_{kl}^* \sum_n\left| \phi_n^X \right\rangle \left\langle\phi_n^X | \theta_i^X \right\rangle \otimes \sum_{\alpha}\left| \psi_{\alpha}^Y \right\rangle \left\langle\psi_{\alpha}^Y | \lambda_j^Y \right\rangle \sum_n \left\langle\theta_k^X | \phi_n^X \right\rangle \left\langle\phi_n^X \right| \otimes \sum_{\alpha} \left\langle\lambda_l^Y | \psi_{\alpha}^Y \right\rangle \left\langle\psi_{\alpha}^Y \right|$$
I know that I can rearrange the terms a little and use the closure relations for the $\theta $ and $\lambda$ basis but these terms dont have the same indicies, therefore I am stuck. I can't obtain equation $(1)$.
How can I do this?