\section*{Supplemental Information}
\paragraph*{Numerical calculations.}
Starting from the master equations, we can formulate equations of motion for the elements of the density matrix 
\begin{equation}
  \varrho_{n,i;m,j}=\langle n,i|\hat \varrho|m,j\rangle,
\end{equation}
with the first index being the cavity photon number and the second being the atomic excitation ($i,j=1,2$). The explicit equations can be written as
\begin{widetext}
\begin{align}
 \dot\varrho_{n,1;m,1}=&-[i\delta_\text c(n-m)+\tfrac{\kappa}{2}(n+m)]\varrho_{n,1;m,1}-ig[\sqrt{n}\varrho_{n-1,2;m,1}-\sqrt{m}\varrho_{n,1;m-1,2}]\nonumber\\
		     &-i\Omega_\text R[\varrho_{n,2;m,1}-\varrho_{n,1;m,2}]+\Gamma\varrho_{n,2;m,2}+\kappa\sqrt{(n+1)(m+1)}\varrho_{n+1,1;m+1,1},\\
\dot\varrho_{n,1;m,2}=&[i(\delta_\text a-(n-m)\delta_\text c)-\tfrac{\Gamma+\kappa(n+m)}{2}]\varrho_{n,1;m,2}-ig[\sqrt{n}\varrho_{n-1,2;m,2} -\sqrt{m+1}\varrho_{n,1;m+1,1}]\nonumber\\
		      &- i\Omega_\text R(\varrho_{n,2;m,2}-\varrho_{n,1;m,1})+\kappa\sqrt{(n+1)(m+1)}\varrho_{n+1,1;m+1,2},\\
  \dot\varrho_{n,2;m,1}=&-[i(\delta_\text a+(n-m)\delta_\text c)+\tfrac{\Gamma+\kappa(n+m)}{2}]\varrho_{n,2;m,1}-ig[\sqrt{n+1}\varrho_{n+1,1;m,1} -\sqrt{m}\varrho_{n,2;m-1,2}]\nonumber\\
		      &- i\Omega_\text R(\varrho_{n,1;m,1}-\varrho_{n,2;m,2})+\kappa\sqrt{(n+1)(m+1)}\varrho_{n+1,2;m+1,1},\\
  \dot\varrho_{n,2;m,2}=&-[i\delta_\text c(n-m)+\Gamma+\tfrac{\kappa}{2}(n+m)]\varrho_{n,2;m,2}-ig[\sqrt{n+1}\varrho_{n+1,1;m,2}-\sqrt{m+1}\varrho_{n,2;m+1,1}]\nonumber\\
		      &-i\Omega_\text R(\varrho_{n,1;m,2}-\varrho_{n,2;m,1})+\kappa\sqrt{(n+1)(m+1)}\varrho_{n+1,2;m+1,2}.
\end{align}
\end{widetext}
We truncate the set of equations at a sufficiently large photon number $N$. By varying $N$, the validity of the calculations can be checked. Using ${\rm Tr} \hat \varrho =1$, we can eliminate one element of the main diagonal, in our case, we chose $\varrho_{0,1;0,1}$.  This introduces an inhomogeneity into the equations, allowing us to calculate the steady state density matrix simply by inverting the matrix of coefficients and multiplying with the inhomogeneity. Finally, the expectation values of interest can be directly obtained.



\paragraph*{Measurement of Squeezing.} The normally ordered variance of a light field is usually measured via  balanced homodyne detection, see e.g.~\cite{WelVo}. In case of single atom fluorescence the complications stem from the small collection efficiency, substantially reducing the effect to be measured. This problem can be resolved by correlation measurements~\cite{vogel,vogel95}, since in this case all correlation functions include the quantum efficiencies only as proportionality factors. One may detect the intensity correlation function
\begin{equation}
  \mathcal G^{(2,2)}(t_1,t_2)=\langle\hat E_1^{(-)}(t_1)\hat E_2^{(-)}(t_2)\hat E_2^{(+)}(t_2)\hat E_1^{(+)}(t_1)\rangle
\end{equation}
by the scheme in Fig.~2 of Ref.~\cite{vogel95}.
Here the indices $1,2$ refer to the different output channels of a beamsplitter. Following~\cite{vogel95}, we consider the difference between correlations at equal times and the corresponding steady state value,
\begin{equation}
  \Delta\mathcal G^{(2,2)}(t)=\mathcal G^{(2,2)}(t,t)-\lim_{\tau\rightarrow\infty}\mathcal G^{(2,2)}(t,t+\tau).
\end{equation}
For (stationary) single atom fluorescence as the signal and a 50/50 beamsplitter we obtain
\begin{equation}
  \Delta\mathcal G^{(2,2)}=-\frac{1}{4}\left(I_{\text{fl}}^2+I_{\text{lo}}\langle:(\Delta\hat E_{\text{fl}})^2:\rangle\right),
\end{equation}
with 'fl' and 'lo' denoting fluorescence and local oscillator, respectively, and $I$ is the intensity. For $I_\text{lo}\gg I_\text{fl}$, the second term is dominant. Squeezing is detected if $\Delta\mathcal G^{(2,2)}>0$.  The effect of our optimization of squeezing, which yields a factor of two, is clearly observed by this method.


\end{document}




