Squeezed vacuum through a beam splitter in the Fock basis Let's say we have a beam splitter like the one shown in the image below, where we have two inputs: a field (generally the field we want to study) in the signal mode $\hat{a}_0$, and a local oscillator field in the probe mode $\hat{a}_1$. The input-output relations of the beam splitter, assuming a certain symmetry, can be described by the parametrization
\begin{equation}
 \begin{pmatrix} \hat{a}_t \\ \hat{a}_r \end{pmatrix} = \begin{pmatrix}
    \cos\theta & e^{i\varphi}\sin\theta \\ -e^{-i\varphi}\sin\theta & \cos\theta \end{pmatrix} \begin{pmatrix} \hat{a}_0 \\ \hat{a}_1 \end{pmatrix},
\tag{1}\label{eq1}
\end{equation}
where, for example, the selection of parameters $\varphi = 0$, $\theta = \pi/2$ would give us a balanced beam splitter of real coefficients. For more information about the math behind beam splitters, one can check the book by Christopher C. Gerrry and Peter L. Knight titled Introductory Quantum Optics, in particular section 6.2 - Quantum mechanics of beam splitters.

Now, let's consider the input state $|i\rangle = |0\rangle_0|z\rangle_1$, where $|z\rangle = \hat{S}(z)|0\rangle = e^{\frac{z^*}{2}\hat{a}^2 - \frac{z}{2}(\hat{a}^{\dagger})^2}|0\rangle$, $z \in \mathbb{C}$, is a squeezed vacuum state. Using the relations from Eq. \eqref{eq1} to write $\hat{a}_1$ and $\hat{a}_1^{\dagger}$ in terms of $\hat{a}_t$, $\hat{a}_t^{\dagger}$, $\hat{a}_r$ and $\hat{a}_r^{\dagger}$, and using the fact that the action of the beam splitter onto $|0\rangle_0|0\rangle_1$ is to simply turn it into $|0\rangle_t|0\rangle_r$, I obtain the output state
\begin{equation}
|f\rangle = e^{\frac{z^*}{2}\left(e^{-2i\varphi}\sin^2\theta \hat{a}_t^2 + \cos^2\theta \hat{a}_r^2 +2e^{-i\varphi}\sin\theta\cos\theta \hat{a}_t\hat{a}_r\right) - H.c.}|0\rangle_t|0\rangle_r,
\tag{2}\label{eq2}
\end{equation}
where $H.c.$ stands for Hermitian conjugate. My question is about how to obtain an expression of the form
\begin{equation}
|f\rangle = \sum_{n,m=0}^{\infty}C_{n,m}|n\rangle_t|m\rangle_r
\tag{3}\label{eq3}
\end{equation}
for the output state. I've noticed that the argument of the exponential in Eq. \eqref{eq2} seems to contain the arguments of both the single-mode squeezing operators for the transmitted and reflected modes ($\hat{S}_i(z) = e^{\frac{z^*}{2}\hat{a}_i^2 - \frac{z}{2}(\hat{a}_i^{\dagger})^2}$, $i=t, r$), as well as the argument for the two-mode squeezing operator acting on both modes ($\hat{S}_{tr}(z) = e^{\frac{z^*}{2}\hat{a}_t\hat{a}_r - \frac{z}{2}\hat{a}_t^{\dagger}\hat{a}_r^{\dagger}}$), and I know how to find the answer for either of those cases (single-mode or two-mode squeezing) separately, as the decomposition in terms of number states for these cases are well known (for instance, see again the book by C. Gerry and P. Knight, sections 7.1 and 7.7).
In fact, I've tried to proceed as it is done in section 7.7 of the book to obtain the expression in terms of number states by solving an eigenvalue problem starting from $\hat{a}_t|0\rangle_t|0\rangle_r$, then applying the Baker-Hausdorff lemma to introduce my exponential operator and finally substituting the decomposition from Eq. \eqref{eq3} to arrive at a recurrence relation for the coefficients. The relation I've found is of the form $k_1\sqrt{n+1}C_{n,n} + k_2\sqrt{n+1}C_{n+1,n+1} + k_3\sqrt{n+2}C_{n+2,n} = 0$, but I don't know how to solve this analytically, or even numerically. Do you know how this could be solved?
If not, is there any other way to get an analytical expression for $|f\rangle$ in terms of number states? Or is there any kind of numerical method I could use to approximate the action of the exponential from Eq. \eqref{eq2} onto the two-mode vacuum?
 A: In case it might be useful to anyone reading this question in the future, I have found the answer in this paper: Generating Schrödinger-cat-like states by means of conditional measurements on a beam splitter. In this paper, they find the output state conditioned to the outcome of a photon number measurement on the output corresponding to the vacuum input. I will explain part of the process, without the actual derivations, but I will give all the necessary ingredients to arrive at the answer.
In the paper, they consider a similar setup to mine, with the addition of a phase for the transmission coefficients $\varphi_T$ that I omitted for reasons of symmetry. Then, they use a group theoretical approach by expressing the transformation from the BS in terms of the generators of the Lie algebra of SU(2) associated with angular momentum, by defining:
\begin{equation} 
\hat{V} = e^{i\varphi\hat{L}_3}e^{-2i\theta\hat{L}_2}e^{-i\varphi\hat{L}_3},
\end{equation}
where
\begin{equation}
\hat{L}_2 = \frac{1}{2i}(\hat{a}_0^{\dagger}\hat{a}_1 - \hat{a}_1^{\dagger}\hat{a}_0), \hspace{1cm} \hat{L}_3 = \frac{1}{2}(\hat{a}_0^{\dagger}\hat{a}_0 - \hat{a}_1^{\dagger}\hat{a}_1),
\end{equation}
which gives the same transformation as the one I present in my question if we do
\begin{equation}
\hat{a}_t = \hat{V}\hat{a_0}\hat{V}^{\dagger}, \hspace{1cm} \hat{a}_r = \hat{V}\hat{a_1}\hat{V}^{\dagger}.
\end{equation}
Note that, in order to have everything in terms of my original variables, I have considered $\varphi_T = 0, \varphi_R = \varphi$, where $\varphi_T$ and $\varphi_R$ are variables used in the paper.
Then, they switch from the Heisenberg picture (which is what we were considering) to the Schrödinger picture, where the transformation is applied to the state instead of the creation and anihilation operators, and the transformation is given by $\hat{\rho}_{out} = \hat{V}^{\dagger}\hat{\rho}_{in}\hat{V}$. They consider the case where $\hat{\rho}_{in} = \hat{\rho}_{in0} \otimes |0\rangle_1\langle0|$ and they derive the general result of applying the beam splitter transformation, and they later condition the state on the fact that $m$ photons are detected in the reflected $r$ output state by doing
\begin{equation}
\hat{\rho}_{out}(m_r) = \frac{\langle m_r | \hat{\rho}_{out}|m_r\rangle}{\mathrm{Tr}_t(\langle m_r | \hat{\rho}_{out}|m_r\rangle)}.
\end{equation}
As an example, they consider the particular case of squeezed vacuum in the first input, $\hat{\rho}_{in} = |z\rangle_0\langle z| \otimes |0\rangle_1\langle0|$, and this is where I got the answer from. I will not replicate all the derivations here, but the final result they obtain is a pure state ($\rho_{out}(m_r) = |\psi_{m_r}\rangle\langle\psi_{m_r}|$) given by:
\begin{equation}
|\psi_{m_r}\rangle = |\psi_{m_r}(\alpha)\rangle = \frac{1}{\sqrt{\mathcal{N}_{m_r}}}\sum_{n_t=0}^{\infty}c_{m_r,n_t}(\alpha)|n_t\rangle,
\end{equation}
where $\alpha$ is a factor that depends on the squezing variables and the BS angles; for $z= re^{i\phi}$, we have that $\alpha = \cos^2(\theta)\tanh(r)e^{i\phi}$. The normalization factor and the complex coefficients $c_{m_r,n_t}$ are given by:
\begin{equation}
c_{m_r,n_t} = \frac{(n_t+m_r)!}{\Gamma[\frac{1}{2}(n_t+m_r)+1]\sqrt{(n_t)!}}\frac{[1+(-1)^{n_t+m_r}]}{2}\left(\frac{1}{2}\alpha\right)^{\frac{n_t+m_r}{2}},
\end{equation}
\begin{equation}
\mathcal{N}_{m_r} = \frac{1}{\sqrt{1-|\alpha|^2}}\left[\frac{|\alpha|^2}{1-|\alpha|^2} \right]^{m_r}\sum_{k=0}^{[m_r/2]}\frac{(m_r!)^2}{(m_r-2k)!(k!)^2(2|\alpha|)^{2k}},
\end{equation}
where $[m_r/2]$ in the summation means the integral part of $m_r/2$. The final ingredient needed is the probability of measuring $m_r$ photons in the reflected output, which comes out to be
\begin{equation}
P(m_r) = \sqrt{\frac{1-|k|^2}{1-|\alpha|^2}}\left[\frac{|\alpha|^2(1-\cos^2(\theta))}{\cos^2(\theta)(1-|\alpha|^2)} \right]^{m_r}\sum_{k=0}^{[m_r/2]}\frac{m_r!}{(m_r-2k)!(k!)^2(2|\alpha|)^{2k}},
\end{equation}
where $k = \tanh(r)e^{i\phi}$.
