|n_1, n_2 \rangle \}$, with $n_1$ and $n_2$ counting the number of photons in the first and second mode, respectively. The modes losses can be treated at the hamiltonian level by considering $$\label{eq:model_tilde}
\tilde{H}= \hat{H} -i\frac{\hbar\gamma_{+}}{2}\hat{a}^{\dagger}_{+}\hat{a}_{+} -i\frac{\hbar\gamma_{-}}{2}\hat{a}^{\dagger}_{-}\hat{a}_{-}\, ,$$ where $\hat{H}$ is given in Eq. (\[eq:model\]). With straightforward substitution, Eq. (\[eq:model\_tilde\]) can be written as $$\tilde{H}=\sum_{i=1,2}[\tilde{\Delta}_{i} \hat{a}_{i}^{\dag}\hat{a}_{i}+
U_{\mathrm{nl}}\hat{a}_{i}^{\dag}\hat{a}_{i}^{\dag}\hat{a}_{i}\hat{a}_{i}]
+ \tilde{J}(\hat{a}_{1}^{\dag}\hat{a}_{2}+\hat a_{2}^{\dag}\hat{a}_{1})+F(\hat{a}_{1}^{\dagger}+\hat{a}_{1}) \, .$$ where we have defined $\tilde{\Delta}_i= \hbar[\Delta\omega_{i} - (\gamma_{+}+\gamma_{-})/4]$ and $\tilde{J}=J- \hbar(\gamma_{+}-\gamma_{-})/4$. With the ansatz $$\begin{aligned}
\vert \psi \rangle=C_{00} \vert 00\rangle &+& C_{10} \vert 10\rangle + C_{01} \vert 01\rangle 
+ C_{20} \vert 20\rangle \nonumber \\
&+& C_{11} \vert 11\rangle + C_{02} \vert 02\rangle + ... \, ,   \end{aligned}$$ the steady state solution can be found by solving the coupled equations for the coefficients $C_{n_1 n_2}$ from $i\hbar\frac{\partial}{\partial t} \vert \psi \rangle  = \tilde{H} \vert \psi\rangle = 0$ $$\begin{aligned}
 & \tilde{\Delta} C_{10}+J C_{01} + F C_{00} = 0 \,  \label{eq:lowF1}\\
 & \tilde{\Delta}C_{01}+ J C_{10}  = 0 \,  \\
 & 2 (\tilde{\Delta} + U_{\mathrm{nl}}) C_{20}+\sqrt{2} J C_{11} + \sqrt{2} F C_{10} = 0 \,  \\
 & 2 \tilde{\Delta} C_{11}+\sqrt{2} J (C_{20}+C_{02})+ F C_{01} = 0 \, , \\
 & 2 (\tilde{\Delta} + U_{\mathrm{nl}}) C_{02}+\sqrt{2} J C_{11}  = 0 \, ,\label{eq:lowF5} \end{aligned}$$ 