I wish to prove the following theorem:
Theorem: For every linear map $\Phi\in T(\mathcal{X},\mathcal{Y})$ there exist sets of operators $\{A_r: r\in\Sigma\}$ and $\{B_r: r\in\Sigma\}$, with $A_r:L(\mathcal{X})\rightarrow L(\mathcal{X},\mathcal{Y})$ and $B_r:L(\mathcal{X})\rightarrow L(\mathcal{X},\mathcal{Y})$ such that \begin{equation} \Phi(X)=\sum_r A_r X B_r^*,~~~~~~\forall X\in L(\mathcal{X}), \end{equation} where $\mathcal{X}$ is a finite-dimensional complex vector space, $L(\mathcal{X})$ is the set of linear operators from $\mathcal{X}\rightarrow\mathcal{X}$ and $T(\mathcal{X},\mathcal{Y})$ is the set of all linear maps from $L(\mathcal{X})\rightarrow L(\mathcal{Y})$.
This is what I have so far as proof:
Let $\{E_{kl}\}$ and $\{E_{mn}\}$ be the standard bases for $L(\mathcal{X})$ and $L(\mathcal{Y})$ respectively. Expanding $X$ in the operator basis for $L(\mathcal{X})$ and using linearity of $\Phi$ we have \begin{equation} \Phi(X)=\sum_{kl}X_{kl}~\Phi(E_{kl}). \end{equation} The above implies that \begin{equation} \Phi(X)=\sum_{kl,mn}X_{kl}~\langle m |\Phi(E_{kl})|n\rangle |m\rangle\langle n|\\ =\sum_{kl,mn}\langle m |\Phi(E_{kl})|n\rangle |m\rangle\langle k|~X~|l\rangle\langle n|. \end{equation} This implies that \begin{equation}\label{eq:vecphi(X)} vec\Big(\Phi(X)\Big)=\sum_{kl,mn}\langle m |\Phi(E_{kl})|n\rangle~vec\Big( |m\rangle\langle k|~X~|l\rangle\langle n|\Big)\notag\\ =\sum_{kl,mn}\langle m |\Phi(E_{kl})|n\rangle |m\rangle\langle k|\otimes|n\rangle\langle l|~vec\big(X\big)\notag\\ =\Phi\Big(vec\big(X\big)\Big). \end{equation} The last equality follows from the linearity of vec and $\Phi$. Thus, \begin{equation} \Phi=\sum_{kl, mn}\phi_{mk,nl}~|m\rangle\langle k|\otimes|n\rangle\langle l|, \end{equation} where \begin{equation} \phi_{mk,nl}=\langle m |\Phi(E_{kl})|n\rangle. \end{equation}
Note that the above equation can also be expressed as \begin{equation} \Phi=\sum_{kl, mn}\phi_{mk, nl}~|mn \rangle\langle kl|\notag\\ =\sum_{kl, mn}\phi_{mk, nl}~vec\Big(E_{mn}\Big)vec\Big(E_{kl}\Big)^*, \end{equation}
As $\Phi\in T(\mathcal{X},\mathcal{Y})$, it admits a singular value decomposition, namely sets of vectors $\{x_1,x_2,...x_r\}\subset\mathcal{X}\otimes\mathcal{X}$ and $\{y_1,y_2,...y_r\}\subset\mathcal{Y}\otimes\mathcal{Y}$ along with a set of positive reals $\{s_1, s_2,...s_r\}$ such that \begin{equation} \Phi=\sum_rs_ry_rx_r^*. \end{equation} The last two equations together imply \begin{equation} \phi_{mk,nl}=\sum_r s_r vec\Big(E_{mn}\Big)^*y_rx_r^*vec\Big(E_{kl}\Big)\notag\\ =\sum_r s_r ~a_{mn}^r ~b_{kl}^{r*}, \end{equation} where $a_{mn}^r=vec\Big(E_{mn}\Big)^*y_r$ and $b_{kl}^{r*}=x_r^*~vec\Big(E_{kl}\Big)$. Recalling the second equation from the top for $\Phi(X)$ and using the expression for $\phi_{mknl}$ we have \begin{equation} \Phi(X)=\sum_{kl,mn}\phi_{mk, nl} |m\rangle\langle k|~X~|l\rangle\langle n|\\ =\sum_r s_r \sum_{kl,mn} ~a_{mn}^r ~b_{kl}^{r*} |m\rangle\langle k|~X~|l\rangle\langle n|. \end{equation}
Now, if one could interchange the indices $n$ and $k$ one would be through as then the sum of $m,k$ would give an operator that would multiply $X$ from the left and summing over $l,n$ would give an operator that would multiply $X$ from the right. However, that is certainly not obvious.
I know an alternative proof that uses the Choi-Jamilkowski isomorphism. However, I'm not sure if this proof can be completed or if it is wrong.
Any comment or reference would be of great help.
Thanks.
I believe that your method of proof will not work. First of all, your Theorem as stated is incorrect, since your operators $A_r$ and $B_r$ are not on the correct spaces. The correctly stated theorem is:
Theorem: For every linear map $\Phi\in \operatorname{T}(\mathcal{X},\mathcal{Y})$ there exist collections of operators $\{A_r: r\in\Sigma\}\subset\operatorname{L}(\mathcal{X},\mathcal{Y})$ and $\{B_r: r\in\Sigma\}\subset\operatorname{L}(\mathcal{X},\mathcal{Y})$ such that $$ \Phi(X)=\sum_r A_r X B_r^*,~~~~~~\forall X\in L(\mathcal{X}). $$
That is, each $A_r$ and $B_r$ is a linear operator from $\mathcal{X}$ to $\mathcal{Y}$ (not from $\operatorname{L}(\mathcal{X})$ to $\operatorname{L}(\mathcal{X},\mathcal{Y})$, as your theorem statement indicates).
It is helpful to introduce the following notation to be able to discriminate between the mapping and its matrix representation (both of which you call $\Phi$ in your question). The natural representation of a linear mapping $\Phi\in\operatorname{T}(\mathcal{X},\mathcal{Y})$ is the operator $K(\Phi)\in\operatorname{L}(\mathcal{X}\otimes\mathcal{X},\mathcal{Y}\otimes\mathcal{Y})$ defined as $$ K(\Phi)x\otimes y = \operatorname{vec}(\Phi(xy^{\scriptscriptstyle\mathsf{T}})) $$ for all vectors $x,y\in\mathcal{X}$ (or analgously as $K(\Phi)\operatorname{vec}(X)=\operatorname{vec}(\Phi(X))$ for all operators $X\in\operatorname{L}(\mathcal{X})$).
The "matrix entries" of the natural representation $K(\Phi)$ are now what you call $\phi_{mk,nl}$. Namely, the values $\phi_{mk,nl}$ satisfy $$ K(\Phi)(\operatorname{vec}(E_{k,l})) = \sum_{m,n}\phi_{mk,nl}\operatorname{vec}(E_{m,n}). $$ Following your attempted proof, we can find a singular value decomposition of $K(\Phi)$ by choosing operators $A_r\in\operatorname{L}(\mathcal{X})$ and $B_r\in\operatorname{L}(\mathcal{Y})$ such that $$ K(\Phi) = \sum_r \operatorname{vec}(B_r)\operatorname{vec}(A_r)^* \tag{$1$} $$ (where we can absorb the values $s_r$ into either $A_r$ or $B_r$), and these are exactly the operators whose matrix entres you have denoted as $a_{kl}^r$ and $b_{mn}^r$ respectively. Unfortunately, this collection of operators $A_r$ and $B_r$ cannot be the desired operators in your Theorem since they are not operators of the correct form. Namely, the operators $A_r$ and $B_r$ in the decomposition in (1) are in $A_r\in\operatorname{L}(\mathcal{X})$ and $B_r\in\operatorname{L}(\mathcal{Y})$. But the operators you are looking for in your theorem are $A_r,B_r\in\operatorname{L}(\mathcal{X},\mathcal{Y})$. This is therefore not the correct singular value decomposition needed to prove the Theorem.
The "correct" way to prove the Theorem is simply to use the Choi-Jamiolkowski representation of the mapping, which is the operator $J(\Phi)\in\operatorname{L}(\mathcal{Y}\otimes\mathcal{X})$ defined as $$ J(\Phi) = \sum_{kl}\Phi(E_{k,l})\otimes E_{k,l}. $$ The desired operators $\{A_r: r\in\Sigma\}\subset\operatorname{L}(\mathcal{X},\mathcal{Y})$ and $\{B_r: r\in\Sigma\}\subset\operatorname{L}(\mathcal{X},\mathcal{Y})$ now come from the singular value decomposition of $J(\Phi)$, where $$ J(\Phi) = \sum_r \operatorname{vec}(A_r)\operatorname{vec}(B_r^{\scriptscriptstyle\mathsf{T}})^*. $$