text:
# Photon Statistics of a Single Atom Intracavity EIT-Kerr System

## I Introduction

A promising avenue of research in the attempt to build quantum computers uses large optical nonlinearities to create essential elements for quantum computation, such as quantum gates Nielsen00 . Obtaining large, noiseless nonlinearities is of course a fundamental goal in the field of nonlinear optics, and technological advances in recent years, for example in cavity quantum electrodynamics (CQED) Berman94  and the ability to access and study strongly coupled quantum systems, offer an exciting new era in this field.

Attempts to achieve large optical nonlinearities are generally plagued by a seemingly necessary trade-off between large absorption and the strength of nonlinearity: to obtain a strong nonlinearity one has to drive the atomic system close to resonance, triggering large absorption and therefore a strong source of noise. However, it is possible to reduce atomic absorption (and hence spontaneous emission) by utilizing quantum coherence effects in multilevel atoms. In particular, electromagnetically induced transparency (EIT) Harris97  was employed by Schmidt and Imamoğlu Schmidt96  to devise a scheme involving four–level atoms which produces a large Kerr nonlinearity with virtually no noise. It was then shown by Imamoğlu et al. Imam97  that if such a strong optical Kerr nonlinearity is implemented in a CQED setting, then it is possible to realize the effect of photon blockade, in which the atom–cavity system effectively acts as a turnstile device for single photons. The physical explanation of this effect is simple. Only the transition between the ground and first excited state of the nonlinear atom–cavity system is resonant with the driving field. Higher states are detuned from resonance by an amount proportional to the strength of nonlinearity. Tian and Carmichael Tian92  have also predicted such an effect with a single two-level atom strongly coupled to the cavity mode.

The proposal of Schmidt and Imamoğlu Schmidt96 , although very appealing in its use of EIT to reduce decoherence, meets obstacles in the attempt for implementation in the many-atom CQED environment. It was shown that this particular realization of the photon blockade system is not ideal for demonstrating photon blockade, because of the strong linear dispersion of the medium Grangier98 ; Gheri99 . In an attempt to avoid this difficulty, Rebić et al. Rebic99  proposed a model in which a single four–level atom is trapped in a high–finesse microcavity. It was shown that this system (called the single-atom EIT-Kerr system) can act as a near-ideal Kerr optical nonlinearity. In such a strongly coupled system, the excitations can be labeled as polaritons, which are defined as mixtures of atom/cavity mode excitations. For weak to moderate driving, the EIT–Kerr system can be approximated by a two-state system, corresponding to the two lowest lying polariton eigenstates. The transitions to the next set of excited states (the second manifold), obtained by introducing a second photon from the driving field into the cavity, are highly detuned from the bare–cavity resonance, and therefore cannot be excited. Hence, in effect, the weakly driven single–atom EIT–Kerr system acts as an ideal photon blockade device.

In a further work Rebic01 , a Hamiltonian for the effective two-level system was derived using a polariton basis, and it was shown that such a Hamiltonian can describe the spectral properties and energy splittings (dynamic Stark effect) seen in the full model. Furthermore, to explain the properties of the second order correlation function, it was shown that the effective two-level system has to be generalized to include two extra states in the first excitation manifold Rebic00 . If more than one atom is introduced the situation becomes more complicated, since each additional atom introduces extra energy levels into the second manifold, some of which are connected to the lower levels via transitions which are resonant with the cavity mode. However, it was shown by Werner and Imamoğlu Werner99  that the introduction of an additional atomic detuning can be used to solve this problem (see also the work of Greentree et al. Green00 ).

In Rebic00 , a comparison of the EIT–Kerr scheme and the extended Jaynes–Cummings scheme (i.e. the standard Jaynes-Cummings model JCmodel  with the effects of losses included) was made in terms of their effectiveness in producing photon blockade. The EIT–Kerr scheme exhibited a value of $`g^{(2)}(0)`$ (the second order correlation function at zero time delay) of $`\frac{1}{40}`$ that of the corresponding Jaynes–Cummings scheme. This was attributed to destructive quantum interference between certain transition amplitudes. In particular, with a suitable choice of parameters, probability amplitudes for transitions from the first excited state (in the first manifold) to the two second manifold eigenstates closest to the cavity resonance cancel each other. This leads to enhanced antibunching of the photons leaving the cavity so that $`g^{(2)}(0)0`$.

We have remarked that the validity of these results depends on the assumption of weak driving. It is therefore of interest to explore how robust the photon blockade is as the driving field strength is increased. We have performed an analysis of the driving field dependence in the EIT–Kerr system, and briefly compared the results with those obtained for the extended Jaynes–Cummings model. A surprising feature is found in the EIT–Kerr system, where a sudden change in photon statistics (i.e. of $`g^{(2)}(0)`$) of light emerging from the cavity occurs at intermediate driving strengths. If plotted against the mean intracavity photon number, $`g^{(2)}(0)`$ is seen to undergo a sudden transition, switching rapidly between antibunching and strong bunching. The exact position of the threshold depends on the characteristic system parameters, namely the atom–field interaction strength and the (EIT) coupling field Rabi frequency.

In Section II we present our model and the methods of solution employed. Section III contains the results of numerical simulations for the photon statistics, showing the surprising behaviour of the second order correlations in the regime of intermediate driving. Section IV contains a detailed explanation for the results presented in the preceding Section, obtained by analyzing the density matrix expressed in an appropriate basis. In Section V we present analytical calculations based on insights obtained in Section IV.

## II Model

### II.1 Hamiltonian

The driven atom-cavity configuration and the atomic energy level scheme are shown schematically in Fig. 1. The Hamiltonian describing the system in the rotating wave and electric dipole approximations is $`=_0+_d`$, where, in the interaction picture,

$`_0`$ $`=`$ $`\mathrm{}\delta \sigma _{22}+\mathrm{}\mathrm{\Delta }\sigma _{44}+i\mathrm{}g_1\left(a^{}\sigma _{12}\sigma _{21}a\right)`$ (1a)
$`+i\mathrm{}\mathrm{\Omega }_c\left(\sigma _{23}\sigma _{32}\right)+i\mathrm{}g_2\left(a^{}\sigma _{34}\sigma _{43}a\right),`$
$`_d`$ $`=`$ $`i\mathrm{}_p\left(aa^{}\right).`$ (1b)
In these equations, $`\sigma _{ij}`$ are atomic raising and lowering operators (for $`ij`$), and atomic energy level population operators (for $`i=j`$); $`a`$ and $`a^{}`$ are the cavity field annihilation and creation operators, respectively. Detunings $`\delta `$ and $`\mathrm{\Delta }`$ are defined from the relevant atomic energy levels; $`g_{1,2}`$ are atom–field coupling constants for the transitions $`|1|2`$ and $`|3|4`$, and $`\mathrm{\Omega }_c`$ is the Rabi frequency of a coherent field coupling the transition $`|2|3`$. The cavity driving field is characterized through the parameter $`_p`$, related to the power output of the driving laser $`𝒫`$ via

$$_p=\sqrt{\frac{𝒫\kappa T^2}{4\mathrm{}\omega _{cav}}},$$
(2)
where $`T`$ is the cavity mirror transmission coefficient, $`\kappa `$ is the cavity decay rate, and $`\omega _{cav}`$ is the cavity mode frequency. Damping due to cavity decay and spontaneous emission is introduced below in the context of the quantum trajectory approach Carmichael93 .

Four atomic levels plus the cavity mode span a Hilbert space of dimension $`4\times N`$, where $`N`$ denotes the truncation order in a Fock state expansion of the cavity field subspace. In the absence of driving, Hamiltonian $`_0`$, given by Eq. (1a), takes a block–diagonal form, with $`N`$ blocks on the main diagonal. Each block represents a manifold of eigenstates associated with the appropriate term in the Fock expansion. The ground, first and second manifolds have been analyzed from the viewpoint of photon blockade in Refs Rebic99 ; Werner99 ; Green00 . The general structure of the dressed states in an arbitrary $`n^{\mathrm{th}}`$ manifold has been discussed in Ref. Rebic01 .

Addition of the driving term (1b) complicates the analysis, since the Hamiltonian matrix does not retain the block–diagonal form. It is possible, however, to re–express the driving Hamiltonian in terms of polariton operators and effective Rabi frequencies for transitions between dressed states. This has been done in Rebic01 , and in the remainder of this article we will draw on these results.

Our analysis of the driven ‘atom-cavity molecule’ will proceed in two complementary directions. First, using the method of quantum trajectories Carmichael93 , we obtain numerical data. Then, using the formalism of Ref. Rebic01 , we construct an effective Hamiltonian in the polariton basis, which encapsulates the physics sufficiently to explain the numerical results.

### II.2 Damping: Quantum Trajectories

In the quantum trajectories approach, damping enters the model through collapse operators, with each of these corresponding to one decay channel. In the case under consideration we require the following four collapse operators,

$`C_1`$ $`=`$ $`\sqrt{\gamma _1}\sigma _{12},C_2=\sqrt{\gamma _2}\sigma _{32},`$
$`C_3`$ $`=`$ $`\sqrt{\gamma _3}\sigma _{34},C_4=\sqrt{\kappa }a,`$ (3)
where $`\gamma _k`$ denote spontaneous emission rates into each of the decay channels, and $`\kappa `$ is the cavity field decay rate. The effective Hamiltonian used in the trajectories approach is non–Hermitian and takes the form

$$_{eff}=i\underset{k=1}{\overset{4}{}}C_k^{}C_k,$$
(4)
with $``$ given by (1).

In deciding on the truncation for the cavity mode Hilbert space, special care must be taken to include a sufficient number of states to capture all the relevant dynamics. If an empty cavity would be driven by an external coherent field corresponding to the parameter $`_p`$, the amplitude of the intracavity coherent field would be

$$\alpha =i_p/\kappa ,$$
(5)
and the expected mean intracavity photon number $`n=|\alpha |^2`$. Our simulations include driving strengths of up to $`_p=3\kappa `$, so the truncation is set at $`N=40`$. The inset in Fig. 2 shows that the actual mean intracavity photon number stays well below its empty cavity counterpart for the given range of driving, thus justifying the chosen truncation.

## III Numerical Simulations of Photon Statistics

In this section we present the results of numerical simulations using the method of quantum trajectories Carmichael93 , with results averaged over 10000 trajectories. Values of parameters chosen for the simulations are consistent with recent experiments Hood00 , so the scheme presented in this paper should be experimentally viable.

We evaluate the second–order correlation function $`g^{(2)}(0)`$ as a function of driving strength $`_p`$. It was established earlier Imam97 ; Rebic99  that this function is a good measure of photon blockade; perfect photon blockade yields perfectly antibunched photons. The steady–state second–order correlation function is given by

$$g^{(2)}(\tau )=\underset{t\mathrm{}}{lim}\frac{a^{}(t)a^{}(t+\tau )a(t+\tau )a(t)}{a^{}(t)a(t)a^{}(t+\tau )a(t+\tau )}.$$
(6)
Perfect antibunching/photon blockade corresponds to $`g^{(2)}(\tau =0)=0`$; for a coherent field $`g^{(2)}(\tau =0)=1`$, and for a bunched field $`g^{(2)}(\tau =0)>1`$ Walls94 .

Fig. 2 shows $`g^{(2)}(\tau =0)`$ as a function of the (scaled) driving parameter for a single–atom EIT–Kerr system, compared with the same quantity calculated for the extended Jaynes–Cummings model. The extended JC model is driven on the lower Rabi resonance, as envisaged by Tian and Carmichael Tian92 . In the weak driving regime, simulation confirms earlier results – stronger photon blockade in the EIT–Kerr system. As the driving increases, the extended Jaynes-Cummings model gradually saturates, and the field correlation asymptotically tends to one. Naively, one would expect qualitatively similar behavior in the EIT–Kerr model. Our simulation, however, shows a vastly different result. Over a narrow range of driving, the statistics of the field changes from strongly antibunched to strongly bunched, and the coherent state value $`g^{(2)}(0)=1`$ is approached asymptotically from above as $`_p`$ is increased further.

Note that the quantity being increased here is the scaled driving parameter, so one may argue that it is the ratio that determines this behavior, i.e. we may either increase the driving strength or decrease the cavity decay rate to obtain the same result. However, to clarify this issue as related to the photon statistics, in Fig. 3 we show a parametric plot of second–order correlation function against the expectation value of intracavity photon number, $`a^{}a_{ss}`$. Both quantities are now unscaled by any external parameter, and the sudden nature of the switch becomes even more obvious, showing a phase–transition–like behaviour. This can be related back to Fig. 2. In particular, the transition from antibunching to bunching happens over a range of driving for which the intracavity photon number stays practically constant (see inset of Fig. 2), leading to the suddenness of the transition seen in Fig. 3, and the concentration of points around the threshold region. To further emphasize this, in Fig. 3 $`(b)`$ we plot the ‘pure’ second–order correlation function $`a^{}a^{}aa_{ss}`$ which also exhibits a threshold-like behaviour.

One other feature of the numerical results is noted. The bunching–antibunching transition is sharper and the subsequent bunching stronger for smaller $`\mathrm{\Omega }_c`$. In fact, it is the increasing ratio $`g_1/\mathrm{\Omega }_c`$ that really matters. The sharper transitions also occur at smaller values of $`a^{}a_{ss}`$. For a decreasing ratio $`g_1/\mathrm{\Omega }_c`$, $`g^{(2)}(0)`$ approaches the behaviour of the extended JC model.

## IV Density Matrix Treatment

We proceed to determine which eigenstates of the strongly coupled quantum system contribute significantly to its dynamics. The total density matrix of the system in the steady state can be written in most general terms as

$$\rho =\underset{a,b}{}c_{ab}|ab|,$$
(7)
where $`a`$, $`b`$ belong to the set of all possible states of the system in an arbitrary basis. The natural basis for the simulation is the one of bare states. The cavity mode subspace is truncated at 40, making the size of $`\rho `$ equal to $`160\times 160`$. Let $`T`$ be the transformation that diagonalizes Hamiltonian (1), i.e. the full Hamiltonian, with driving included. The density matrix can be transformed into a new basis as $`\rho _{diag}=T\rho T^1`$, and we can look for the nonzero elements of this matrix. Diagonal elements of the matrix correspond to populations of the dressed states, while off-diagonal elements correspond to coherences between the dressed states. The results depend on the size of $`_p`$, i.e. we expect the number of non–zero elements to increase as $`_p`$ is increased. For a large part of the $`_p`$ range considered, however, the density matrix is dominated by elements corresponding to two particular states.

Figure 4 shows the nonzero matrix elements of density matrix $`\rho _{diag}`$. Diagonal states are sorted in increasing order; the state with largest negative detuning is at $`(i,j)=(1,\mathrm{\hspace{0.17em}1})`$, while the state with largest positive detuning is at $`(i,j)=(160,\mathrm{\hspace{0.17em}160})`$. States with the smallest detuning (i.e. closest to the cavity resonance) are at the centre of matrix, at positions 80 and 81 along the main diagonal. Two diagonal elements dominate the matrix, and we identify these as being the populations of the Stark-split states $`|\psi _\pm `$ (see Appendix A and Rebic01 ). Stark-split states are therefore found at the positions 80 and 81 on the main diagonal of the density matrix. At positions 79 and 82 are the two states from the second manifold (two-photon excitations), closest to the resonance. At the positions 78 and 83 are two remaining states from the first manifold; at the positions 77 and 84 are two states from the third manifold (three-photon excitations), closest to the resonance, and so on. The states we have just identified suffice to indicate the dynamics of the system. Off-diagonal elements are coherences between the appropriate dressed states.

One striking feature can be noted immediately from Fig. 4. Namely, the first ‘square’ of elements (4 elements on each side) encircling the centre square (Stark states, $`2\times 2`$) remains much smaller than all of the other accessible states. In terms of dressed states, this means that the second manifold states remain unpopulated and the coherences between these and all the other states vanish as well. Furthermore, as the driving increases, the states whose coherences with the Stark states increase to a noticeable size are the third manifold states (see Fig. 4 for $`_p=0.09\kappa `$). As the driving increases further, the remaining two first manifold states (their populations and coherences with the Stark states) also start to contribute. So, contrary to expectations, the contribution to the dynamics of the states beyond the effective two-level system do not increase according to the smallness of their detunings from the bare cavity frequency. The two least–detuned second-manifold states are in fact essentially unpopulated.

From the above analysis of the density matrix elements, we can deduce the solution to the ‘photon statistics puzzle’ of Figs. 2 and 3. The strong photon antibunching at small $`_p`$ is the consequence of the effective two-level system, and has been well understood. The sharp rise in $`g^{(2)}(0)`$ with increasing $`_p`$ can be attributed to the two-photon transitions needed to populate the third manifold states from the first manifold, and then decay back in cascade to the first manifold. These two-photon decays cause the sharp increase in the second order output field correlations. But, what causes the system to skip second manifold states? Again, the answer can be deduced from the density matrix. Strong coherence between the first and third manifold states, together with the missing population in the second manifold states (and vanishing coherences associated with these states) uncovers the effect of EIT-type quantum interference between the dressed states. This is not surprising, since the quantum interference between the transitions from first to second manifold have been already discussed in Ref. Rebic00 .

These features are shown in Fig. 5, where the relevant populations and coherences are shown. Note that the combined population of the two inner second manifold states (Fig. 5 $`(b)`$, solid line) vanishes for a large interval of driving, since the value of $`10^{16}`$ is at the numerical precision value, and fluctuations are numerical, not physical in nature. These populations become nonzero at the value of driving strength at which $`g^{(2)}(0)`$ of Fig. 2 (solid line) peaks. These plots further justify the discussion of the preceeding paragraphs.

We proceed with the development of an effective model with relatively few levels which nevertheless captures most of the important features of the dynamics.

## V Effective Model

In the formulation of an effective model, we rely on the formalism developed in Ref. Rebic01 . This formalism was very successful in the development of an effective two-level theory, and in explaining the fluorescence spectrum. Now we extend the model and include a total of six dressed states in the effective model. These states are shown in Fig. 6.

It should be noted, however, that the dressed states shown in Fig. 6 do not correspond exactly to the dressed states discussed in Section IV. Namely, the dressed states of Sec. IV are often referred to as doubly dressed states, since they diagonalize the Hamiltonian with driving contributions included. In this Section, we will treat driving separately (see Rebic01 ), and let the dressed states represent the eigenstates of the interaction Hamiltonian (1a) alone. Driving can then be included through effective Rabi frequencies $`\mathrm{\Omega }_{ij}`$, coupling dressed states $`|e_i`$ and $`|e_j`$.

The relevant effective non-Hermitian Hamiltonian (including driving and damping) is thus

$`_{eff}`$ $`=`$ $`\mathrm{}ϵ_2p_2^{}p_2+\mathrm{}ϵ_3p_3^{}p_3+\mathrm{}ϵ_4p_{j4}^{}p_{j4}+\mathrm{}ϵ_5p_{j5}^{}p_{j5}`$ (8)
$`+i\mathrm{}\mathrm{\Omega }_{01}\left(p_1p_1^{}\right)`$
$`+i\mathrm{}\left(\mathrm{\Omega }_{12}^{}p_2\mathrm{\Omega }_{12}p_2^{}\right)+i\mathrm{}\left(\mathrm{\Omega }_{13}^{}p_3\mathrm{\Omega }_{13}p_3^{}\right)`$
$`+i\mathrm{}\mathrm{\Omega }_{24}\left(p_{24}p_{24}^{}\right)+i\mathrm{}\mathrm{\Omega }_{25}\left(p_{25}p_{25}^{}\right)`$
$`+i\mathrm{}\mathrm{\Omega }_{34}\left(p_{34}p_{34}^{}\right)+i\mathrm{}\mathrm{\Omega }_{35}\left(p_{35}p_{35}^{}\right)`$
$`i\mathrm{}\mathrm{\Gamma }_0p_1^{}p_1i\mathrm{}\mathrm{\Gamma }_{22}p_2^{}p_2i\mathrm{}\mathrm{\Gamma }_{33}p_3^{}p_3`$
$`i\mathrm{}\mathrm{\Gamma }_{44}p_{j4}^{}p_{j4}i\mathrm{}\mathrm{\Gamma }_{55}p_{j5}^{}p_{j5}`$
$`i\mathrm{}\mathrm{\Gamma }_{23}p_2^{}p_3i\mathrm{}\mathrm{\Gamma }_{32}p_3^{}p_2`$
$`i\mathrm{}\mathrm{\Gamma }_{45}p_{j4}^{}p_{j5}i\mathrm{}\mathrm{\Gamma }_{54}p_{j5}^{}p_{j4}.`$
Operators $`p_{jk}`$ are polariton operators defined by $`|e_k=p_{jk}^{}|e_j`$. For the states $`|e_k`$ which are accessible from only one lower state, the notation has been abbreviated, so, for example, $`p_2p_{12}`$. Note also that the indices $`j`$ are dummy indices, i.e. $`p_{jk}^{}p_{jl}=|e_ke_l|`$. Rabi frequencies $`\mathrm{\Omega }_{12}`$ and $`\mathrm{\Omega }_{13}`$ have a phase term making them purely imaginary, while all the other Rabi frequencies are real. Damping terms were discussed in detail in Rebic01 . It is now possible to formalize the distinction between the two dressed states bases used in this Section and Section IV. The basis in which the Hamiltonian (8) is written is defined in Ref. Rebic01  and depicted in Fig. 6. The discussion of Section IV is based upon diagonalizing the total Hamiltonian (including all of the 160 states used for the numerical simulation).

With the effective Hamiltonian (8), we can also write the master equation for an effective density matrix as

$`\dot{\rho }_{eff}`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\left(_{eff}\rho _{eff}\rho _{eff}_{eff}^{}\right)`$ (9)
$`+2{\displaystyle \underset{i,j}{}}S_i\rho _{eff}S_j^{},`$
where $`S_k`$ now denote the polariton collapse operators (Appendix B.1). The effective density matrix $`\rho _{eff}`$ has the dimension $`6\times 6`$ \- a significant reduction from $`160\times 160`$ used to obtain the results in Sec. IV. Equations of motion for the elements of $`\rho _{eff}`$ are given in Appendix B.2.

From the equations of motion, we can uncover terms which lead to the effect of quantum interference. For example, equations for the populations of the second manifold states $`\rho _{22}`$ and $`\rho _{33}`$, depend on the populations $`\rho _{44}`$ and $`\rho _{55}`$ of the third manifold states and the coherences between these states $`\rho _{45}`$ and $`\rho _{54}`$; the latter with rate $`\mathrm{\Gamma }_{45}`$. At the same time, equations for $`\rho _{44}`$ and $`\rho _{55}`$ do not depend on the second manifold states, nor their mutual coherence. The same holds for the coherences $`\rho _{23}`$ and $`\rho _{45}`$ and their adjoints; the equation for $`\rho _{23}`$ depends on $`\rho _{44}`$, $`\rho _{55}`$, $`\rho _{45}`$ and $`\rho _{54}`$, but not vice versa. Population of and coherence between the second manifold states is linked to the population of and coherence between the third manifold states. If this dependence is removed from the equations of motion, i.e., terms dependent on $`\rho _{44}`$, $`\rho _{55}`$, $`\rho _{45}`$ and $`\rho _{54}`$ are removed from the equations for $`\dot{\rho }_{22}`$, $`\dot{\rho }_{33}`$, $`\dot{\rho }_{23}`$ and $`\dot{\rho }_{32}`$, cancellation of the population in the second manifold ceases to occur.

The 35 equations of the effective model can, in principle, be solved analytically in the steady state. However, the resulting expressions are complicated and do not offer significant physical insight, so we have opted to perform numerical solutions of the equations given in Appendix B.2 and check for validity of the effective model. Having the solutions for the populations and coherences, the second order correlation function for zero time delay can be calculated as the ratio of

$`a^{}a^{}aa=|w_{01}|^2|w_{12}|^2\rho _{22}+|w_{01}|^2|w_{13}|^2\rho _{33}`$ (10a)
$`+[|w_{12}|^2|w_{24}|^2+|w_{13}|^2|w_{34}|^2+(w_{12}^{}w_{24}^{}w_{13}w_{34}+\mathrm{c}.\mathrm{c}.)]\rho _{44}`$
$`+[|w_{12}|^2|w_{25}|^2+|w_{13}|^2|w_{35}|^2+(w_{13}^{}w_{35}^{}w_{12}w_{25}+\mathrm{c}.\mathrm{c}.)]\rho _{55}`$
$`+(|w_{01}|^2w_{12}w_{13}\rho _{34}+\mathrm{c}.\mathrm{c}.)`$
$`+[(|w_{12}|^2w_{24}w_{25}^{}+|w_{13}|^2w_{34}w_{35}^{}`$
$`+w_{12}^{}w_{25}^{}w_{13}w_{34}+w_{12}w_{24}w_{13}^{}w_{35}^{})\rho _{45}+\mathrm{c}.\mathrm{c}.],`$
and the square of
$`a^{}a=|w_{01}|^2\rho _{11}+|w_{12}|^2\rho _{22}+|w_{13}|^2\rho _{33}`$ (10b)
$`+\left(|w_{24}|^2+|w_{34}|^2\right)\rho _{44}+\left(|w_{25}|^2+|w_{35}|^2\right)\rho _{55}`$
$`+[w_{12}w_{13}^{}\rho _{23}+(w_{24}w_{25}^{}+w_{34}w_{35}^{})\rho _{45}+\mathrm{c}.\mathrm{c}.]`$
where $`w_{ij}=\mathrm{\Omega }_{ij}/_p`$, and c.c. stands for complex conjugate.

The results are displayed in Fig. 7. We have compared the numerical solutions of the effective model with the results of Sec. III and found a very good qualitative agreement. We do find strong bunching and threshold behaviour occurring at the same values of $`_p`$ and $`a^{}a_{ss}`$. However, the effective model differs from the full simulations in the size of $`g^{(2)}(0)`$ at its peak, by approximately a factor of 2. Including more states in the effective model would yield full agreement with the numerical data. In particular, we found that including two states closest to the resonance from up to eight manifolds reproduces the numerical data exactly. The reason is that the two-photon cascade decay responsible for the behaviour of $`g^{(2)}(0)`$ can result from the decay of the higher states to third manifold states first, i.e. two-photon cascade can, in the manner of speaking, be driven “from below” and “from above”. Naturally, the decay of higher lying states introduces more single-photon transitions as well. Therefore, adding one manifold at a time to an effective model reveals that the increase in maximum value of $`g^{(2)}(0)`$ gradually diminishes with new manifolds added, settling at its maximum value after the inclusion of the eighth manifold. On the other end, the effective model of Fig. 6 is the smallest possible model which (at least qualitatively) reproduces the strong bunching effect in this system.

## VI Dynamics of the Forward Scattering of Light

To obtain a different and useful perspective on the physical processes involved in the changing nature of the statistics of light emitted by the coupled atom-cavity system, we can split the field operator into contributions from a coherent mean amplitude and from an incoherent part Carmichael85 ,

$$a=\alpha +\mathrm{\Delta }a,$$
(11a)
where $`\alpha a`$ denotes the coherent amplitude component of the intracavity field, while $`\mathrm{\Delta }a`$ denotes the incoherent amplitude component $`\mathrm{\Delta }a=0`$, whose emergence is the result of scattering of the cavity field by the atom. Using this decomposition, one can rewrite the expression for $`g^{(2)}(0)`$ as
$`g^{(2)}(0)1`$ $`=`$ $`{\displaystyle \frac{a^{}a^{}aa_{ss}}{a^{}a_{ss}^2}}`$ (11b)
$`=`$ $`{\displaystyle \frac{:(\alpha ^{}\mathrm{\Delta }a+\alpha \mathrm{\Delta }a^{})^2:_{ss}}{\left(|\alpha |^2+\mathrm{\Delta }a^{}\mathrm{\Delta }a_{ss}\right)^2}}`$
$`+{\displaystyle \frac{4|\alpha |\mathrm{Re}\left[\mathrm{\Delta }a^2\mathrm{\Delta }a_{ss}\right]}{\left(|\alpha |^2+\mathrm{\Delta }a^{}\mathrm{\Delta }a_{ss}\right)^2}}`$
$`+{\displaystyle \frac{\mathrm{\Delta }a^2\mathrm{\Delta }a^2_{ss}\left(\mathrm{\Delta }a^{}\mathrm{\Delta }a_{ss}\right)^2}{\left(|\alpha |^2+\mathrm{\Delta }a^{}\mathrm{\Delta }a_{ss}\right)^2}}`$
$`=`$ $`S(\mathrm{\Delta }a)+T(\mathrm{\Delta }a)+V(\mathrm{\Delta }a),`$
where $`::`$ denotes normal ordering, and $`a^{}a_{ss}=|\alpha |^2+\mathrm{\Delta }a^{}\mathrm{\Delta }a_{ss}`$. The three terms in this expansion, denoted by S, T and $`V`$, have been identified by Carmichael Carmichael85  for the case of a two-level atom. The decomposition (11) shows how the behaviour of $`g^{(2)}(0)`$ for different values of driving field can be interpreted as the effect of self-homodyning between the coherent and incoherent components of the intracavity field Mandel82 . From this viewpoint, it is easy to identify $`S(\mathrm{\Delta }a)`$ as a term describing the squeezing in the field quadrature in phase with the driving field, $`V(\mathrm{\Delta }a)`$ gives the variance in the incoherent component, and $`T(\mathrm{\Delta }a)`$ describes intensity-amplitude correlations in the incoherent component. Both $`V(\mathrm{\Delta }a)`$ and $`T(\mathrm{\Delta }a)`$ are determined by the correlations in the intensity, so the departure from the coherent value of the correlation function can be assigned to the effects of squeezing and the effects of intensity correlations.

The contributions from the incoherent component of the intracavity field are shown in Fig. 8. The squeezing and intensity correlation parts are shown separately. The antibunching for weak fields comes from the squeezed fluctuations, which reduce the contribution from the coherent scattering. As the driving increases, the squeezing decreases, but the variance in intensity fluctuations becomes negative, so the remaining antibunching comes from the sub-Poissonian intensity fluctuations in the incoherent component of the field. For $`_p\kappa `$, the squeezing contribution effectively vanishes, while $`T`$ and $`V`$ become positive, and antibunching disappears. Strong bunching clearly originates in the super-Poissonian intensity correlations, and the correlation function is dominated by $`V`$.

## VII Nonclassical Behaviour of the Correlation Function

The presence of nonclassical effects in a driven atom-cavity system is a topic that has received much attention for many years. The system usually studied has been the canonical system of quantum optics – a single two-level atom confined in an optical cavity. Photon statistics in the bad-cavity limit was thoroughly studied by Rice and Carmichael Rice88 , who analyzed the sub-Poissonian statistics and photon antibunching in the cavity transmission, for the case of weak driving. Their analysis was extended by Carmichael et al. Carmichael91  to a system containing $`N`$ two-level atoms. This analysis was further refined by Brecha et al. Brecha99 . Clemens and Rice Clemens00  have extended the consideration involving a single atom to include arbitrary driving field strength and dephasing. In their analysis, Clemens and Rice pay special attention to nonclassical effects known as ‘undershoots’ and ‘overshoots’. These are related to the violation of inequalities that hold for classical correlations, in particular violations that occur not in the value of $`g^{(2)}(\tau =0)`$, but for certain time delays $`\tau >0`$. The explanation for the undershoots has been given by Carmichael et al. Carmichael91  in terms of quantum interference of probability amplitudes and collapse of the wavefunction.

The requirements for the classicality of the field correlations can be derived from the Cauchy-Schwartz inequality (see Brecha99  and references therein), and expressed in terms of the second-order correlation function as

$$|g^{(2)}(\tau )1||g^{(2)}(0)1|.$$
(12)
Values in excess of those allowed classically are called overshoots, while values below are called undershoots. Overshoots have been observed recently by Mielke et al. Mielke98 .

In this context, it is of interest to see if the overshoots and/or undershoots can be found in the single-atom EIT-Kerr system under consideration. Photon antibunching, as an example of nonclassical photon statistics, has already been predicted Rebic99 , and the effective two-level behaviour analyzed Rebic99 ; Werner99 ; Rebic01 . In the present article, we have shown how the effects of self-homodyning of squeezed dipole radiation yields photon antibunching in the low to moderate driving limit. We have also shown that quantum interference between the probability amplitudes contributes to both strong antibunching and strong bunching, for weak and strong driving fields, respectively. Given this range of behaviours, we might therefore expect undershoots and overshoots to also occur in the single-atom EIT-Kerr system under suitable conditions.

Fig. 9 shows correlation functions for several values of driving field strength. The values of driving have been chosen where nonclassical behaviour is expected to be found. For weak driving, where the antibunching is strong, the delay-time dependence of the correlation function is well-understood in terms of the effective two-level system. The interesting region is for those values of driving for which $`g^{(2)}(0)`$ increases through 1, the value for a coherent field. This is also the region in which the dynamics is well-described in terms of quantum interference and increased incidence of two-photon emissions.

Not surprisingly, this is also the range of parameters where largest violations of the classical inequalities occur. Since the increase in the value of $`g^{(2)}(0)`$ is due to the purely quantum effect of interference between probability amplitudes, non-classical behaviour of the correlations can be expected. The explanation of these effects given by Rice and Carmichael Rice88  and Carmichael et al. Carmichael91 , although in a different context, still holds. As shown in Sec. VI, self-homodyning of squeezed dipole radiation with the driving field occurs in the EIT-Kerr system in a similar manner to that for a two-level atom. An alternative explanation in terms of quantum interference of the driving field with the atomic polarization after the collapse of the wave function upon a photon detection event offers even more insight. This is best understood in the context of quantum trajectory theory.

This explanation is illustrated in Fig. 10, where single trajectory realizations for the intracavity field and photon number are plotted. We see that the collapses tend to occur in succession before the system returns to a (quasi-) steady state (which for this region of parameters happens after time $`\gamma _j^1`$). Once it returns into a steady state, a few Rabi cycles pass before the next set of collapses occur. The value of $`g^{(2)}(0)`$ is determined by the ratio of the number of jumps upwards to number of jumps downwards in photon number, where jump upwards suggests that the detection of a photon increases the probability of detecting a second photon immediately afterwards. Naturally, at the value of driving where $`g^{(2)}(0)`$ peaks (see Fig. 2), collapses are almost exclusively upwards, as illustrated in Figure 11. Undershoots appear as the consequence of a change in sign that the amplitude undergoes at the collapse Carmichael91 . As the system returns towards its steady state the polarization becomes close in magnitude and opposite in sign to the driving field, producing a near-zero mean intracavity field, leading to the reduced detection probability for a second photon.

The overshoots can also be explained in terms of the collapses of the wave function. The detection of the first photon, emitted from the steady state situation, collapses the wave function of the system. The subsequent time evolution as the system returns to the steady state determines the photon correlations. For weak driving field, the probability for the second collapse to occur before the system returns to the steady state is extremely small, since it is proportional to the mean intracavity photon number. For stronger driving fields the probability for subsequent collapses increases, specifically due to the large correlations between first and third manifold states, as shown in Fig. 5. Therefore, the probability for a second photon detection after some time $`\tau <\gamma ^1`$ increases as well, causing the correlation function overshoot. In experiment, such event pairs are the source of delayed coincidence counts. A third collapse is also likely to occur before the steady state is reached. However, the overshoot disappears (or significantly decreases) for all time delays $`\tau `$ after the first peak. This is expected, since the exact form of the wave function after the second collapse depends on the delay time between the second and third photon. Averaging over all possible realizations washes out the nonclassical effects due to different possible evolutions following the second collapse.

A stronger driving field causes more subsequent collapses to occur, and nonclassical correlations are completely washed out. We find that overshoots and undershoots vanish at driving strength $`_p1.14\kappa `$.

## VIII Conclusion

We have presented an analysis of the properties of the photon statistics of light emitted by a single atom intracavity EIT-Kerr system. It was found that the statistics change abruptly as the driving field strength increases. Specifically, strong photon antibunching, dominant in the weak to moderate driving regimes, is replaced by a strong photon bunching in the output field for the stronger driving. We have identified the effect of quantum interference between the dressed states to be responsible for this sudden change, and presented an effective model explaining qualitative features of this behaviour.

Furthermore, we have analyzed contributions from the incoherent scattering to the system dynamics and found that the strong photon antibunching can be explained in terms of the self-homodyning of the incoherent intracavity component with the coherent component for weak driving, and in terms of reduced intensity fluctuations for moderate driving strengths. Strong bunching is the signature of super-Poissonian intensity fluctuations.

Finally, nonclassical behaviour of the correlation function was found, and the effects of undershoots and overshoots analyzed.

###### Acknowledgements.

The authors would like to thank M. J. Collett and M. K. Olsen for many valuable comments and suggestions. This work was supported by the Marsden Fund of the Royal Society of New Zealand and The University of Auckland Research Committee.

## Appendix A Stark Splitting of the Dressed States

In this Appendix, we briefly review the effect of dynamic Stark splitting of dressed states, fully elaborated upon in Ref. Rebic01 .

Dressed states of the Hamiltonian (1a) contain, among others, two states on exact cavity resonance, separated by energy $`\mathrm{}\omega _{cav}`$. These are the ground state and a state belonging to the first excitation manifold. In the schematic depiction of Fig. 6, these are denoted by $`|e_0`$ and $`|e_1`$. Driving by the external field $`_p`$ enters through the effective Rabi frequency $`\mathrm{\Omega }_{01}=_p/\sqrt{1+g_1^2/\mathrm{\Omega }_c^2}`$. Resonant coupling of the two levels causes splitting of the two states into a symmetric and antisymmetric linear combination of the two

$$|\psi _\pm =\left(|e_0\pm |e_1\right)/\sqrt{2},$$
(13)
with corresponding energies

$$ϵ_\pm =\pm \sqrt{\mathrm{\Omega }_{01}^2\left(\mathrm{\Gamma }_0/2\right)^2},$$
(14)
where $`\mathrm{\Gamma }_0=\kappa /\left(1+g_1^2/\mathrm{\Omega }_c^2\right)`$ is the decay rate of the excited state $`|e_1`$. The splitting of energy levels occurs at the threshold values of driving field

$$_p=\frac{\kappa /2}{\sqrt{1+g_1^2/\mathrm{\Omega }_c^2}}.$$
(15)
The two states $`|\psi _\pm `$ are the two states closest to the cavity resonance for the range of driving strengths considered in this article.

## Appendix B Quantum Jumps and Density Matrix Elements

In this Appendix, we derive jump terms for the effective master equation (9), and the equations of motion for density matrix elements.

### B.1 Jump Terms

We follow the notation of Carmichael Carmichael93  and rewrite the effective master equation (9) in a Lindblad form as $`\dot{\rho }_{eff}=_{eff}\rho _{eff}`$, where the Liouvillian superoperator can be divided into a part describing the free evolution between the jumps, $`_{eff}𝒮_{eff}`$, and a part describing the jumps, $`𝒮_{eff}`$. Jump term $`𝒮_{eff}\rho _{eff}`$ arises from the equivalent jump term in the full master equation,

$`𝒮\rho `$ $`=`$ $`2\kappa a\rho a^{}+2\gamma _1\sigma _{12}\rho \sigma _{21}`$ (16)
$`+2\gamma _2\sigma _{32}\rho \sigma _{23}+2\gamma _3\sigma _{34}\rho \sigma _{43}.`$
Using the methods of Ref. Rebic01 , and consistent with the truncation of the dressed states space (as in Fig. 6), we can write field and atomic collapse operators in terms of the bare states as

$`a`$ $`=`$ $`|0,\mathrm{\hspace{0.17em}1}1,\mathrm{\hspace{0.17em}1}|+|0,\mathrm{\hspace{0.17em}3}1,\mathrm{\hspace{0.17em}3}|+|0,\mathrm{\hspace{0.17em}4}1,\mathrm{\hspace{0.17em}4}|`$ (17a)
$`+\sqrt{2}\left(|1,\mathrm{\hspace{0.17em}1}2,\mathrm{\hspace{0.17em}1}|+|1,\mathrm{\hspace{0.17em}2}2,\mathrm{\hspace{0.17em}2}|+|1,\mathrm{\hspace{0.17em}3}2,\mathrm{\hspace{0.17em}3}|\right)`$
$`+\sqrt{3}|2,\mathrm{\hspace{0.17em}1}3,\mathrm{\hspace{0.17em}1}|,`$
$`\sigma _{12}`$ $`=`$ $`|1,\mathrm{\hspace{0.17em}1}1,\mathrm{\hspace{0.17em}2}|+|2,\mathrm{\hspace{0.17em}1}2,\mathrm{\hspace{0.17em}2}|,`$ (17b)
$`\sigma _{32}`$ $`=`$ $`|1,\mathrm{\hspace{0.17em}3}1,\mathrm{\hspace{0.17em}2}|+|2,\mathrm{\hspace{0.17em}3}2,\mathrm{\hspace{0.17em}2}|,`$ (17c)
$`\sigma _{34}`$ $`=`$ $`|0,\mathrm{\hspace{0.17em}3}0,\mathrm{\hspace{0.17em}4}|+|1,\mathrm{\hspace{0.17em}3}1,\mathrm{\hspace{0.17em}4}|.`$ (17d)
It can be deduced from the bare states that the collapse operators can be expressed as linear combinations of the following polariton operators: $`a\{p_1,p_2,p_3,p_{24},p_{34},p_{25},p_{35}\}`$; $`\sigma _{12},\sigma _{34}\{p_2,p_3,p_{24},p_{34},p_{25},p_{35}\}`$ and $`\sigma _{32}\{p_2^{}p_2,p_3^{}p_3,p_4^{}p_4,p_5^{}p_5,p_2^{}p_3,p_3^{}p_2,p_{j4}^{}p_{j5},p_{j5}^{}p_{j4}\}`$. Note that the polariton jump terms arising from the atomic jumps associated with the operator $`\sigma _{32}`$ (i.e. proportional to $`\gamma _2`$) couple dressed states within the same manifold. This is possible since the associated atomic transition is not coupled to the cavity mode, but to the classical field $`\mathrm{\Omega }_c`$, and therefore jumps occurring in this atomic decay channel do not switch between the adjacent manifolds.

Using the correspondences above, one can write all of the jump terms in the polariton basis, thus making the transition $`𝒮\rho 𝒮_{eff}\rho _{eff}`$. This transformation can be viewed as a rotation of a truncated basis of a superoperator. The coefficients $`\mathrm{\Gamma }_{ijkl}`$ of the jump terms are given in Table 1.

### B.2 Equations of Motion

The equations of motion for the density matrix elements in the basis schematically shown in Fig. 6 can be found from the master equation (9). There is a total of 35 equations, since population conservation can be used to eliminate one of them. The equations for populations are

$`\dot{\rho }_{00}`$ $`=`$ $`2\mathrm{\Gamma }_0\rho _{11}+\mathrm{\Omega }_0\left(\rho _{01}+\rho _{10}\right),`$ (18a)
$`\dot{\rho }_{11}`$ $`=`$ $`2\mathrm{\Gamma }_0\rho _{11}+\mathrm{\Gamma }_{1122}\rho _{22}+\mathrm{\Gamma }_{1133}\rho _{33}`$ (18b)
$`\mathrm{\Omega }_0\left(\rho _{01}+\rho _{10}\right)`$
$`+\mathrm{\Omega }_{12}\rho _{12}+\mathrm{\Omega }_{12}^{}\rho _{21}+\mathrm{\Omega }_{13}\rho _{13}+\mathrm{\Omega }_{13}^{}\rho _{31},`$
$`\dot{\rho }_{22}`$ $`=`$ $`2\left(\mathrm{\Gamma }_{22}+\mathrm{\Gamma }_{2222}\right)\rho _{22}`$ (18c)
$`\mathrm{\Gamma }_{23}\left(\rho _{32}+\rho _{23}\right)+2\mathrm{\Gamma }_{45}\left(\rho _{45}+\rho _{54}\right)`$
$`\mathrm{\Omega }_{12}\rho _{12}\mathrm{\Omega }_{12}^{}\rho _{21}+\mathrm{\Gamma }_{2233}\rho _{33}+\mathrm{\Gamma }_{2244}\rho _{44}+\mathrm{\Gamma }_{2255}\rho _{55}`$
$`+\mathrm{\Omega }_{24}\left(\rho _{24}+\rho _{42}\right)+\mathrm{\Omega }_{25}\left(\rho _{25}+\rho _{52}\right),`$
$`\dot{\rho }_{33}`$ $`=`$ $`2\left(\mathrm{\Gamma }_{33}+\mathrm{\Gamma }_{3333}\right)\rho _{33}`$ (18d)
$`\mathrm{\Gamma }_{23}\left(\rho _{32}+\rho _{23}\right)+2\mathrm{\Gamma }_{45}\left(\rho _{45}+\rho _{54}\right)`$
$`\mathrm{\Omega }_{13}\rho _{13}\mathrm{\Omega }_{13}^{}\rho _{31}+\mathrm{\Gamma }_{3322}\rho _{22}+\mathrm{\Gamma }_{3344}\rho _{44}+\mathrm{\Gamma }_{3355}\rho _{55}`$
$`+\mathrm{\Omega }_{34}\left(\rho _{34}+\rho _{43}\right)+\mathrm{\Omega }_{35}\left(\rho _{35}+\rho _{53}\right),`$
$`\dot{\rho }_{44}`$ $`=`$ $`2\left(\mathrm{\Gamma }_{44}+\mathrm{\Gamma }_{4444}\right)\rho _{44}+\mathrm{\Gamma }_{4455}\rho _{55}2\mathrm{\Gamma }_{45}\left(\rho _{45}+\rho _{54}\right)`$ (18e)
$`\mathrm{\Omega }_{24}\left(\rho _{24}+\rho _{42}\right)\mathrm{\Omega }_{34}\left(\rho _{34}+\rho _{43}\right),`$
$`\dot{\rho }_{55}`$ $`=`$ $`2\left(\mathrm{\Gamma }_{55}+\mathrm{\Gamma }_{5555}\right)\rho _{55}+\mathrm{\Gamma }_{5544}\rho _{44}2\mathrm{\Gamma }_{45}\left(\rho _{45}+\rho _{54}\right)`$ (18f)
$`\mathrm{\Omega }_{25}\left(\rho _{25}+\rho _{52}\right)\mathrm{\Omega }_{35}\left(\rho _{35}+\rho _{53}\right).`$
The equations for coherences are

$`\dot{\rho }_{01}`$ $`=`$ $`\mathrm{\Gamma }_0\rho _{01}+2\mathrm{\Gamma }_{02}\rho _{12}+2\mathrm{\Gamma }_{03}\rho _{13}`$ (19a)
$`+\mathrm{\Omega }_0\left(\rho _{11}\rho _{00}\right)+\mathrm{\Omega }_{12}\rho _{02}+\mathrm{\Omega }_{13}\rho _{03}`$
$`\dot{\rho }_{02}`$ $`=`$ $`\left(\mathrm{\Gamma }_{22}+\mathrm{\Gamma }_{2222}iϵ_2\right)\rho _{02}\mathrm{\Gamma }_{23}\rho _{03}`$ (19b)
$`+\mathrm{\Gamma }_{0214}\rho _{14}+\mathrm{\Gamma }_{0215}\rho _{15}`$
$`+\mathrm{\Omega }_0\rho _{12}\mathrm{\Omega }_{12}^{}\rho _{01}+\mathrm{\Omega }_{24}\rho _{04}+\mathrm{\Omega }_{25}\rho _{05}`$
$`\dot{\rho }_{03}`$ $`=`$ $`\left(\mathrm{\Gamma }_{33}+\mathrm{\Gamma }_{3333}iϵ_3\right)\rho _{03}\mathrm{\Gamma }_{23}\rho _{02}`$ (19c)
$`+\mathrm{\Gamma }_{0314}\rho _{14}+\mathrm{\Gamma }_{0315}\rho _{15}`$
$`+\mathrm{\Omega }_0\rho _{13}\mathrm{\Omega }_{13}^{}\rho _{01}+\mathrm{\Omega }_{34}\rho _{04}+\mathrm{\Omega }_{35}\rho _{05}`$
$`\dot{\rho }_{04}`$ $`=`$ $`\left(\mathrm{\Gamma }_{44}+\mathrm{\Gamma }_{4444}+iϵ_4\right)\rho _{04}\mathrm{\Gamma }_{45}\rho _{05}`$ (19d)
$`\mathrm{\Omega }_{24}\rho _{02}\mathrm{\Omega }_{34}\rho _{03}+\mathrm{\Omega }_0\rho _{14}`$
$`\dot{\rho }_{05}`$ $`=`$ $`\left(\mathrm{\Gamma }_{55}+\mathrm{\Gamma }_{5555}+iϵ_5\right)\rho _{05}\mathrm{\Gamma }_{45}\rho _{04}`$ (19e)
$`\mathrm{\Omega }_{25}\rho _{02}\mathrm{\Omega }_{35}\rho _{03}+\mathrm{\Omega }_0\rho _{12}`$
$`\dot{\rho }_{12}`$ $`=`$ $`\left(\mathrm{\Gamma }_0+\mathrm{\Gamma }_{22}+\mathrm{\Gamma }_{2222}iϵ_2\right)\rho _{12}\mathrm{\Gamma }_{23}\rho _{13}`$ (19f)
$`+\mathrm{\Gamma }_{1224}\rho _{24}+\mathrm{\Gamma }_{1225}\rho _{25}+\mathrm{\Gamma }_{1234}\rho _{34}+\mathrm{\Gamma }_{1235}\rho _{35}`$
$`+\mathrm{\Omega }_{12}^{}\left(\rho _{22}\rho _{11}\right)`$
$`\mathrm{\Omega }_0\rho _{02}+\mathrm{\Omega }_{13}^{}\rho _{32}+\mathrm{\Omega }_{24}\rho _{14}+\mathrm{\Omega }_{25}\rho _{15}`$
$`\dot{\rho }_{13}`$ $`=`$ $`\left(\mathrm{\Gamma }_0+\mathrm{\Gamma }_{33}+\mathrm{\Gamma }_{3333}iϵ_3\right)\rho _{13}\mathrm{\Gamma }_{23}\rho _{12}`$ (19g)
$`+\mathrm{\Gamma }_{1324}\rho _{24}+\mathrm{\Gamma }_{1325}\rho _{25}+\mathrm{\Gamma }_{1334}\rho _{34}+\mathrm{\Gamma }_{1335}\rho _{35}`$
$`+\mathrm{\Omega }_{13}^{}\left(\rho _{33}\rho _{11}\right)`$
$`\mathrm{\Omega }_0\rho _{03}+\mathrm{\Omega }_{12}^{}\rho _{23}+\mathrm{\Omega }_{34}\rho _{14}+\mathrm{\Omega }_{35}\rho _{15}`$
$`\dot{\rho }_{14}`$ $`=`$ $`\left(\mathrm{\Gamma }_0+\mathrm{\Gamma }_{44}+\mathrm{\Gamma }_{4444}iϵ_4\right)\rho _{14}\mathrm{\Gamma }_{45}\rho _{15}\mathrm{\Omega }_0\rho _{04}`$ (19h)
$`+\mathrm{\Omega }_{24}\rho _{12}+\mathrm{\Omega }_{34}\rho _{13}+\mathrm{\Omega }_{12}^{}\rho _{24}+\mathrm{\Omega }_{13}^{}\rho _{34}`$
$`\dot{\rho }_{15}`$ $`=`$ $`\left(\mathrm{\Gamma }_0+\mathrm{\Gamma }_{55}+\mathrm{\Gamma }_{5555}iϵ_5\right)\rho _{15}\mathrm{\Gamma }_{45}\rho _{14}\mathrm{\Omega }_0\rho _{05}`$ (19i)
$`+\mathrm{\Omega }_{25}\rho _{12}+\mathrm{\Omega }_{35}\rho _{13}+\mathrm{\Omega }_{12}^{}\rho _{25}+\mathrm{\Omega }_{13}^{}\rho _{35}`$
$`\dot{\rho }_{23}`$ $`=`$ $`\left[\mathrm{\Gamma }_{22}+\mathrm{\Gamma }_{2222}+\mathrm{\Gamma }_{33}+\mathrm{\Gamma }_{3333}+\mathrm{\Gamma }_{2233}i\left(ϵ_2+ϵ_3\right)\right]\rho _{23}`$ (19j)
$`\mathrm{\Gamma }_{23}\left(\rho _{22}+\rho _{33}\right)+\mathrm{\Gamma }_{2344}\rho _{44}+\mathrm{\Gamma }_{2355}\rho _{55}+\mathrm{\Gamma }_{2332}\rho _{32}`$
$`+\mathrm{\Gamma }_{2345}\rho _{45}+\mathrm{\Gamma }_{2354}\rho _{54}+\mathrm{\Omega }_{12}^{}\rho _{13}+\mathrm{\Omega }_{13}^{}\rho _{21}`$
$`+\mathrm{\Omega }_{34}\rho _{24}+\mathrm{\Omega }_{35}\rho _{25}+\mathrm{\Omega }_{24}\rho _{43}+\mathrm{\Omega }_{25}\rho _{53}`$
$`\dot{\rho }_{24}`$ $`=`$ $`\left[\mathrm{\Gamma }_{22}+\mathrm{\Gamma }_{2222}+\mathrm{\Gamma }_{44}+\mathrm{\Gamma }_{4444}+\mathrm{\Gamma }_{2244}i\left(ϵ_2+ϵ_4\right)\right]\rho _{24}`$ (19k)
$`\mathrm{\Gamma }_{23}\rho _{34}\mathrm{\Gamma }_{45}\rho _{25}+\mathrm{\Gamma }_{2435}\rho _{35}`$
$`\mathrm{\Omega }_{24}\left(\rho _{44}\rho _{22}\right)\mathrm{\Omega }_{12}\rho _{14}\mathrm{\Omega }_{34}\rho _{23}+\mathrm{\Omega }_{25}\rho _{54}`$
$`\dot{\rho }_{25}`$ $`=`$ $`\left[\mathrm{\Gamma }_{22}+\mathrm{\Gamma }_{2222}+\mathrm{\Gamma }_{55}+\mathrm{\Gamma }_{5555}+\mathrm{\Gamma }_{2255}i\left(ϵ_2+ϵ_5\right)\right]\rho _{25}`$ (19l)
$`\mathrm{\Gamma }_{23}\rho _{35}\mathrm{\Gamma }_{45}\rho _{24}+\mathrm{\Gamma }_{2534}\rho _{34}`$
$`\mathrm{\Omega }_{25}\left(\rho _{55}\rho _{22}\right)\mathrm{\Omega }_{12}\rho _{15}\mathrm{\Omega }_{35}\rho _{23}+\mathrm{\Omega }_{24}\rho _{45}`$
$`\dot{\rho }_{34}`$ $`=`$ $`\left[\mathrm{\Gamma }_{33}+\mathrm{\Gamma }_{3333}+\mathrm{\Gamma }_{44}+\mathrm{\Gamma }_{4444}+\mathrm{\Gamma }_{3344}i\left(ϵ_3+ϵ_4\right)\right]\rho _{34}`$ (19m)
$`\mathrm{\Gamma }_{23}\rho _{24}\mathrm{\Gamma }_{45}\rho _{35}+\mathrm{\Gamma }_{3425}\rho _{25}`$
$`\mathrm{\Omega }_{34}\left(\rho _{44}\rho _{33}\right)\mathrm{\Omega }_{13}\rho _{14}+\mathrm{\Omega }_{35}\rho _{54}\mathrm{\Omega }_{24}\rho _{32}`$
$`\dot{\rho }_{35}`$ $`=`$ $`\left[\mathrm{\Gamma }_{33}+\mathrm{\Gamma }_{3333}+\mathrm{\Gamma }_{55}+\mathrm{\Gamma }_{5555}+\mathrm{\Gamma }_{3355}i\left(ϵ_3+ϵ_5\right)\right]\rho _{35}`$ (19n)
$`\mathrm{\Gamma }_{23}\rho _{25}\mathrm{\Gamma }_{45}\rho _{34}+\mathrm{\Gamma }_{3524}\rho _{24}`$
$`\mathrm{\Omega }_{35}\left(\rho _{55}\rho _{33}\right)\mathrm{\Omega }_{13}\rho _{15}+\mathrm{\Omega }_{34}\rho _{45}\mathrm{\Omega }_{25}\rho _{32}`$
$`\dot{\rho }_{45}`$ $`=`$ $`\left[\mathrm{\Gamma }_{44}+\mathrm{\Gamma }_{4444}+\mathrm{\Gamma }_{55}+\mathrm{\Gamma }_{5555}+\mathrm{\Gamma }_{4455}i\left(ϵ_4+ϵ_5\right)\right]\rho _{45}`$ (19o)
$`\mathrm{\Gamma }_{45}\left(\rho _{44}+\rho _{55}\right)+\mathrm{\Gamma }_{4554}\rho _{55}`$
$`\mathrm{\Omega }_{24}\rho _{25}\mathrm{\Omega }_{34}\rho _{35}\mathrm{\Omega }_{25}\rho _{42}\mathrm{\Omega }_{35}\rho _{43}.`$
Damping coefficients $`\mathrm{\Gamma }_0`$ and $`\mathrm{\Gamma }_{ij}`$ have been evaluated in general in Ref. Rebic01 . Coefficients $`\mathrm{\Gamma }_{ijkl}`$ are related to the jump operators and are given in Table 1. The source of the different terms in these equations should be apparent from the earlier discussion. For example, damping emerging from jump terms is denoted as $`\mathrm{\Gamma }_{ijkl}`$, while damping coming from the effective Hamiltonian (8) is denoted as $`\mathrm{\Gamma }_0`$ and $`\mathrm{\Gamma }_{ij}`$.

It is straightforward to prove that the following equalities hold:

$`2\left(\mathrm{\Gamma }_{22}+\mathrm{\Gamma }_{2222}\right)`$ $`=`$ $`\mathrm{\Gamma }_{1122}+\mathrm{\Gamma }_{3322},`$ (20a)
$`2\left(\mathrm{\Gamma }_{33}+\mathrm{\Gamma }_{3333}\right)`$ $`=`$ $`\mathrm{\Gamma }_{1133}+\mathrm{\Gamma }_{2233},`$ (20b)
$`2\left(\mathrm{\Gamma }_{44}+\mathrm{\Gamma }_{4444}\right)`$ $`=`$ $`\mathrm{\Gamma }_{2244}+\mathrm{\Gamma }_{3344}+\mathrm{\Gamma }_{5544},`$ (20c)
$`2\left(\mathrm{\Gamma }_{55}+\mathrm{\Gamma }_{5555}\right)`$ $`=`$ $`\mathrm{\Gamma }_{2255}+\mathrm{\Gamma }_{3355}+\mathrm{\Gamma }_{4455},`$ (20d)
preserving the population conservation condition $`_j\rho _{jj}=1`$, i.e. $`_j\dot{\rho }_{jj}=0`$.