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}$.
Best Answer
The input coherent state is : $$|\alpha\rangle\otimes|\beta\rangle = \exp\left(-\frac{|\alpha|^2+|\beta|}{2}\right)\exp\left(\alpha \hat a_i^\dagger + \beta \hat b_i^\dagger\right)|0\rangle$$
Then, since $\hat T$ is unitary and $\hat T^\dagger|0\rangle=|0 \rangle $ we have :
\begin{align} \hat T|\alpha\rangle\otimes|\beta\rangle &= \exp\left(-\frac{|\alpha|^2+|\beta|}{2}\right) \exp\left(\hat T\left(\alpha \hat a_i^\dagger + \beta \hat b_i^\dagger\right)\hat T^\dagger\right)|0\rangle \\ &=\exp\left(-\frac{|\alpha|^2+|\beta|}{2}\right) \exp\left(\alpha \hat a_o^\dagger + \beta \hat b_o^\dagger\right)|0\rangle \\ &=\exp\left(-\frac{|\alpha|^2+|\beta|}{2}\right) \exp\left( (\alpha\cos(\theta)-i\beta\sin(\theta))\hat a_i^\dagger + (\beta\cos(\theta) - i\alpha\sin(\theta)) \hat b_i^\dagger\right)|0\rangle \\ &=|\alpha\cos(\theta)-i\beta\sin(\theta)\rangle\otimes|\beta\cos(\theta) - i\alpha\sin(\theta) \rangle \end{align}