Monte-Carlo Simulations of Superradiant Lasing

We present simulations of the superradiant dynamics of ensembles of atoms in the presence of collective and individual atomic decay processes. We unravel the density matrix with Monte-Carlo wave-functions and identify the quantum jumps in a reduced Dicke state basis, which reflects the permutation symmetry of the identical atoms. While the number of density matrix elements in the Dicke representation increases polynomially with atom number, the quantum jump dynamics populates only a single Dicke state at the time and thus efficient simulations can be carried out for tens of thousands of atoms. The calculated superradiance pulses from initially excited atoms agree quantitatively with recent experimental results with strontium atoms but rapid atom loss in these experiments does not permit steady-state superradiance. By introducing an incident flux of new atoms, the system can maintain a large average atom number, and our theoretical calculations predict lasing with millihertz linewidth despite rapid atom number fluctuations.

Introduction Superradiance is caused by the collective interaction of atoms with a common radiation field and has been the subject of continuous interest since the early proposal by Dicke [1], (see also [2,3]). Recent experiments with tens of thousands of atoms in optical cavities [4][5][6] thus explore the possibility of achieving superradiant lasing in a bad cavity. Opposite to the Schawlow-Townes limit of a normal laser, it has a linewidth set by the rate of energy transfer from single atoms to the cavity mode and can in principle reach millihertz level [7]. However, so far, such performance has been hindered by fast atom loss depleting the atomic ensemble before the coherence is established [4]. To solve this problem, we might feed new atoms to the system to stabilize the atom number. This will invariably introduce noise in the system but the coherence shared among the atoms might still be preserved over the time scale where all atoms have been replaced.
In this Letter, we perform numerical studies of the system depicted in Fig. 1 (a), where atoms interact collectively with a lossy cavity mode while being subject to individual decay, dephasing and excitation (pumping) due to the interaction with their local environment, as well as to loss and feeding. Theoretically, this system can be studied with either the laser master equation or, by elimination of the bad cavity, an atomic superradiance master equation [8]. However, since the number of density matrix elements scales exponentially with the atom number, the calculations in the atomic product states basis are normally restricted to only tens of atoms [9]. Assuming identical interactions and permutation symmetry among all the atoms, calculations can be carried out in a Dicke states basis [10][11][12] or collective numbers basis [13][14][15] (exploiting SU(4) group theory [16]), where the complexity scales only cubically with the atom number and thus permits simulations for hundreds of atoms. In this Letter we employ the Monte Carlo wave-function method (MCWF) to unravel the master equation in the Dicke states basis. This, on the one hand, provides novel insights into the excitation dynamics of the system and, on the other hand, permits simulations for tens of thousands * yzhang@phys.au.dk † iyxz@phys.au.dk ‡ moelmer@phys.au.dk  Figure 1. System and processes. Panel (a) shows two-level atoms located within an optical cavity in the presence of the loss of atoms from the cavity mode volume and the feeding of new atoms into the volume. Panels (b-g) show Dicke state diagrams for four atoms illustrating quantum jumps due to collective decay (b); decay, pumping, and dephasing of individual atoms (c-e); atom loss and feeding (f,g). The thickness of the arrows indicates the probabilities of the different jumps as discussed in the text. of atoms as encountered in experiments.
The MCWF method [17,18] utilizes ensembles of wavefunctions instead of a density matrix to represent a quantum system, and applies random quantum jumps to describe the effect of dissipation. The jumps can also describe back-action associated with detection, and thus form an essential component of so-called quantum trajectories [19] with applications in quantum measurement and control [20]. When a system is allowed to explore all states of a high dimensional Hilbert space, the conditional wave-functions are lower dimensional objects than density matrices. This fact renders the MCWF method an efficient numerical tool for complex system simulations. Furthermore, the symmetries mentioned above restrict the dynamics to sub-spaces and thus can potentially speed up the wave-function calculations.
The theoretical treatment of N atoms that start in the same pure state and experience only collective decay, is restricted to the sub-space of fully symmetric Dicke states |J, M with collective spin quantum numbers J = N/2 and M = −J, ...J representing the total number of excited atoms J + M . The individual decay, dephasing and incoherent excitation as well as the atom loss and feeding seemingly break the symme-try among the atoms. However, if all atoms undergo these processes with the same rates, the density matrix retains its symmetry under permutation of the atoms, but will explore the Dicke states |J, M with J = N/2, N/2 − 1, ...1/2 or 0, and M = −J, ...J [21]. In the following, we review briefly the master equation for the density matrix in the Dicke states basis and then demonstrate its unraveling and effective simulation with the MCWF method.
Superradiance Master Equation and Monte Carlo Wave-Function Simulations The master equation for the density matrix of N two-level atoms with identical transition frequency ω a can be written as where the Hamiltonian H a = (ω a /2) σ z is characterized by the collective Pauli operator σ z = N k=1 σ z k . The cavity mode with frequency ω c and loss rate γ c provides a channel for the collective decay, and in the bad cavity limit the elimination of the cavity field operators yields [8] the collective Lamb shift and the collective decay rate, respectively, and are determined by the detuning χ = ω c − ω a , the cavity loss γ c and the atom-cavity mode coupling g. The last term in Eq. (1) represents the dissipation of individual atoms by three Lindblad terms ρ with a decay rate γ l , a pumping rate κ l and a dephasing rate d l . To realize the incoherent pumping with κ l and atomic population inversion, the atoms may be excited from the lower level via a higher excited level from which they decay rapidly to the upper level of our effective two-level systems [7,9]. Further below, we shall generalize our analysis to the case with varying atom number.
Dicke states |JM are common eigenstates of the collective operators j 2 = l=x,y,z j l 2 and j z [1] with eigenvalues J(J + 1) and M , where j x = (σ + + σ − )/2, j y = i(σ − − σ + )/2 and j z = σ z /2. For J < N/2, the states with same J and M are degenerate, and may be labeled by an additional quantum number α = 1, ..., d J N , which specifies d J N = N ! (2J + 1) / [(N/2 − J)! (N/2 + J + 1)!] different symmetric linear combinations of atomic product states [22]. Due to the identical dissipation of all the atoms, it was shown in [10,11] that while the dissipative processes of a single atom break the symmetry of the many body state, the sum of these processes over all the atoms in Eq. (1) preserves the symmetry and populates states with different α evenly. These processes do not build coherence between states of different α or J, and the equality of density matrix elements ρ JM α,JM α for different α motivates us to introduce a single term ρ J M M = (1/d J N ) α ρ JM α,JM α to represent all of these elements. The resulting effective density matrix is normalized as Tr(ρ) ≡ J,M d J N ρ J M M = 1. Notice that the degeneracy factor appears also in the evaluation of physical observables.
The equations for all the specified density matrix elements are derived in [10,11]. Now we proceed to unravel them with the MCWF method. To this end we introduce an ensemble of wave-functions Here, the effective states |JM are not real states, but symbolic constructions, which disregard the degeneracy of the subspace with same J and M , but allow sampling of the density matrix element ρ J M M as 1 The normalization of the wave-functions is J,M d J N |C i JM | 2 = 1. The absence of coherence between different J allows us to consider only the matrix elements and thus amplitudes of the wave-functions for each separate value of J.
In the Appendices, we detail the evaluation of the Monte Carlo wave-functions. Here, we emphasize the main procedures and highlight the related quantum jumps. The wavefunctions are represented as state vectors in the Dicke states space and are propagated with a non-Hermitian Hamiltonian including the atomic Hamiltonian and the anti-commutator of the dissipative superoperators. The C i JM amplitudes follow separate equations including a reduction of norm associated with different physical processes that we simulate by discrete quantum jumps. The collective decay applies the jump operator σ − on the state-vectors, which leads to the jumps to the Dicke states with the same J but reduced M , cf. Fig. 1 For the individual decay, we use the Clebsch-Gordan (CG) expansion [10,11] to write the Dicke states of N atoms as the product of those of N − 1 atoms and a single atom, and then apply the operator σ − k on the single atom state. Finally, we use the inverse CG expansion to convert the state back to the Dicke states basis of all N atoms. This introduces jumps to Dicke states with M reduced by unity and J changed by s = 0, +1, −1. The branching ratio of the jumps are given in the Table I in the Appendices and favors jumps to states with reduced values of J, cf. Fig. 1(c). We apply the same treatment to the individual pumping and dephasing by applying the operators σ + k and σ z k to the Dicke states, yielding jumps to states with quantum numbers J + s and M + 1 and M , cf. Fig. 1 (d,e) and the Table I in the Appendices. We now consider the quantum jump description of loss of atoms and feeding of new atoms into the system. For that purpose, we assume the atom loss and feeding rates, γ loss and κ f eed , are independent of the atomic internal state. The average atom number follows the rate equation ∂N /∂t = −γ loss N + κ f eed , which yields a steady state Poisson distribution of atom number with mean N s = κ f eed /γ loss . To simulate the loss of an atom we use the CG expansion to decouple the Dicke states of N atoms as a product of states of N − 1 atoms and of a single atom, and then we discard the state of the single atom by a partial trace, simulated by quantum jumps from Dicke states of N atoms to those of N − 1 atoms with probabilities proportional to squared CG coefficients, cf. from the Dicke state of N atoms to one of N +1 atoms, cf. Fig.  1 (g) and the Appendices. After the implementation of any of the above quantum jumps, the state vectors are re-normalized, and the wave-function propagation and further random jumps proceed.
The number of the effective Dicke states of N atoms is (N + 3) (N + 1) /4 for odd N and (N + 2) 2 /4 for even N . However, if the system has no initial coherence between states of different J, such a coherence will never appear since all the jumps happen either within the same ladder or from one to a neighboring ladder, cf. Fig.1 (b-g). Thus, our wave-functions can at any time be expanded on only a single ladder of states, restricting the dimension to 2J +1 ≤ N +1 elements. If the system is initially prepared in a single Dicke state, e.g. the ground or fully-excited state, all coherences between states with different M vanish for all times, and the system evolution can be effectively simulated as an incoherent jump process between states |JM . In this case, the dynamics of a stochastic wavefunction can be visualized as a trajectory in a (N, J, M ) coordinate system, and we only need to evaluate the different jump probabilities to simulate the random evolution of N (t), J(t), and M (t). By simulation of multiple such trajectories we can evaluate ensemble average of physical observables, such as the population of the Dicke states P i JM = |C i JM (t)| 2 , and the radiation through the cavity mode (defined with the collective operators σ + σ − ), and the radiation into the free space

Superradiant Dynamics of a Fixed Number of Atoms
We first illustrate our simulations with a small system of fifty atoms (in the absence of atom loss and feeding), which are initially all excited, i.e. we occupy the Dicke state |J = 25, M = 25 , cf. Fig. 2 (a,b), and initially all in the ground state, i.e. we occupy the Dicke state |J = 25, M = −25 , but pumped incoherently, cf. Fig. 2 (c,d). The main panels show the time evolution of the emitted radiation intensity, while the insets show the quantum jump trajectories among Dicke states.
To study the influence of individual decay with rate γ l , we set γ l = Γ c in Fig. 2 (a) and γ l = 10Γ c in Fig. 2 (b). For weak individual decay, we observe that the system remains close to the fully symmetric Dicke states with J = N/2, while for larger individual decay, it departs from these states and explores states with lower values of J. In all cases, all atoms eventually end up in the ground state |J = 25, M = −25 in the lower right corner in both insets. This dynamics is also reflected in the radiative emission by the system. Figure 2 (a) shows a minor radiation into the free-space (the blue curve), and a dominant radiation through the cavity mode, i.e. superradiance. The different gray curves reflect shot-to-shot variations in the superradiance from the system, which are rather large due to the small number of atoms in our simulations. The average superradiance (the black solid curve) agrees well with the master equation results (the red dashed curve), verifying our MCWF method. Figure 2 (b) shows that for strong individual decay, the superradiance is suppressed and weaker than the radiation into the free-space (the blue curve), which follows the usual exponential decay law.
To study the effect of individual pumping with rate κ l , we set κ l = Γ c in Fig. 2(c), and κ l = 10Γ c in Fig. 2(d). For weak pumping the trajectories start from the lower right corner and propagate along the lower boundary of the Dicke ladders, and finally wander around the left corner with small values of J, M , while the radiation through the cavity mode fluctuates around a steady-state mean value 25Γ c . For strong pumping the trajectories explore more states and end up in a region around J = 15 and M = 8. As a consequence of the higher symmetry of the states and their higher degree of excitation, the collectively enhanced emission probabilities are higher and both the fluctuations and the averaged radiation increase. For even stronger pumping, we find that the trajectories propagate to the upper right corner of the inset Dicke state diagram, and the radiation through the cavity mode is reduced. In all cases, the averaged radiation trajectories are in perfect agreement with the exact results calculated with the master equation. The reduced radiation at strong pumping can be ascribed to the effect of the pump-induced noise [7], and may also be understood as a suppression of the superradiance transition matrix elements ∝ A JM − = (J + M )(J − M + 1) when the atoms do not occupy both the ground and excited states.
Superradiant Lasing After having verified our method with a toy model calculation for tens of atoms, we now turn to the simulation of much larger ensembles, cf. Fig. 3 and 4, as studied by Norcia et. al. [4]. In their experiment, more than 10 5 strontium atoms are trapped in an optical lattice inside an optical cavity and coupled to a cavity mode through their ultra-narrow S 1 0 − P 3 0 optical clock transition. For this problem, the many independent trajectories explore ranges of J and M according to the same diffusion-like process as observed in Fig. 2. However, with N > 10 5 we have verified that the relative fluctuations among trajectories remain very small, and hence we can simulate the system and obtain the information of interest with only few hundred quantum trajectories. In Fig. 3 we show the superradiance pulses for different numbers of initial atoms N (0) (from 1 × 10 5 to 2 × 10 5 ) in the presence of atom loss. The inset shows that the atom number decreases exponentially (the blue dotted curve), which reduces the population of the upper (lower) level before (after) the superradiance pulse. The calculated pulses agree well with the experimental results (the noisy curves) with a fitted atom loss rate of γ loss = 5.58 Hz (compatible with the magnitude estimated in [4]) and negligible dephasing rate d l . If we increase d l (up to γ loss ), the calculated pulses are reduced and shifted to earlier time and thus do not match the experimental results. If we increase γ l by ten times, the calculated pulses are not affected so much. These results of our analysis confirm the assessment in [4] that the system dynamics is dominated by the collective decay and the atom loss process.
The above simulation indicates that all the atoms are lost in about 0.2 second and thus there is not time for these atoms to establish the coherence leading to the lasing with millihertz line-width. However, by feeding new atoms, we can achieve steady-state superradiance. To calculate the corresponding spectrum, we apply the quantum regression theorem and calculate the spectrum with the Fourier transform of the two-time correlation function of the collective atomic raising and lowering operators, S (ω) ∝ Γ c Re´∞ 0 dτ e −iωτ σ + (τ ) σ − (0) . We follow [17] to calculate the correlation function with the following procedure: we first propagate the stochastic wavefunctions to long time (to approach steady-state), and then apply four different combinations of the collective lowering operator and the unit operator on the final wave-functions to initialize four ancillary wave-functions that we propagate with the MCWF method to obtain correlation function trajectories.  Integrated intensity (blue triangles, right axis) and linewidth (red squares, left axis) of the steady-state spectrum from Ns = 10 5 strontium atoms versus pumping rate κ l without (left) and with (right) atom loss and feeding (of atoms in the lower state). κ f eed = γ loss Ns and other parameters are the same as in Fig. 4.
The average of many such trajectories yields the exact correlation function. For more details, please refer to the Appendices. Note that during the evolution the ancillary wave-functions occupy a coherent superposition of only two states (with identical J and M differing by unity). They can thus be propagated by solution of only two coupled equations, and their quantum jumps are also readily implemented and retain their simple two-components form at all times.
The left panel of Fig. 4 shows the intensity and line-width of the steady-state superradiant spectrum from 10 5 strontium atoms, which are pumped individually with increasing rate κ l , in the absence of atom loss and feeding. We see that the intensity increases linearly with the pumping, while the linewidth decreases and approaches a minimum and then increases due to the pump-induced noise. The intensity agrees qualitatively with the prediction in [7], while the linewidth broadening occurs for much smaller pumping because of the negligible dephasing in our simulations. Note that the minimal linewidth achieved is about 16 mHz and may be much smaller for systems with larger atom number.
The right panel of Fig. 4 shows the results in the presence of atom loss and feeding (in the lower state) being able to maintain the large steady-state atom number, N s = 10 5 with κ f eed = γ loss N s . As each atom typically stays for only a fraction of a second in the cavity, the incoherent pumping has to be fast to ensure excited state population available for the lasing process. Our calculations, indeed, show a threshold effect, requiring a pumping rate κ l > 2 Hz to yield a collective emission signal. Since the adjacent Dicke states in the ancillary wavefunctions are exposed to the same jumps during the correlation function evolution, their coherence survives longer than the atoms, leading to sub-hertz linewidth for pumping smaller than 9 Hz and a minimal 60 mHz line-width for the 2 Hz pumping rate.
Conclusions We have developed a Monte Carlo wavefunction (MCWF) method to study superradiance of two-level atoms in a cavity. Our treatment incorporates the atomic collective decay, individual decay, pumping and dephasing as quantum jumps within or among Dicke ladders of fixed atom number, and the atom loss and feeding as jumps between Dicke states with reduced and increased atom number. The wavefunction populates only one Dicke state at a time and the above processes together determine how the wave-function explores the Dicke ladders. Our method is verified by the comparison with the exact master equation calculations for tens of atoms and our calculated superradiance pulses by more than 10 5 strontium atoms agree with the experimental results with a fitted atom loss rate. Atom loss prevents steady-state lasing, but by feeding ground-state atoms into the cavity, we can stabilize the mean atom number and obtain steady-state spectra with a minimal linewidth about 60 mHz. Notice that the corresponding coherence time is about 16 second and is thus hundred times longer than the time spent by any typical atom in the cavity. The MCWF unraveling of the density matrix as illustrated here can provide not only numerical results for superradiant lasing (in cross-over regime [23]), but may also be applied to study superradiant beats [24], synchronization [25], and spin-squeezing [26] of atomic ensembles. Furthermore, by visualizing the dynamical evolution of the quantum states, it can provide important insight into the underlying physical processes.
respectively. The factors A JM ± = (J ∓ M ) (J ± M + 1) in Eq. (A1) depend on the Dicke state quantum numbers. These numbers govern the so-called no-jump dynamics of the wave function. The imaginary part reduces the norm of the state vector by a sum of infinitesimal probabilities β δp β , which are, in the master equation, compensated by (sandwich) feeding terms, e.g., In the Monte Carlo wave-function method, the feeding terms are represented by the return of the population by quantum jumps into different final states. β enumerates eighteen different quantum jump channels, depicted in Fig. 1 (b-g) in the main text and with the probabilities summarized in Table I. The values of the different δp β and the corresponding quantum jump actions, will be detailed in the following. For the collective decay, we apply the collective jump operator σ − acting on the wave-function to yield |ψ (t + δt) = √ δtΓ c σ − |ψ (t) , which leads to the relation for the state amplitudes This quantum jump, illustrated in Fig. 1 (b) in the main text, which lowers M and retains J, occurs with the probability δp c = δtΓ c d J N M |A JM − C JM | 2 . Since the system has vanishing coherences between states with different J, we can restrict the expansion of the quantum state to a single value of J, which is unchanged by the collective decay process (and by the no-jump dynamics).
The decay of individual atoms lowers M by unity and changes J by 0 or ±1 as sketched in Fig. 1 (c) in the main text. It is represented by the master equation terms . Their effect on the states, k σ q k |JM JM | (σ q k ) + , is evaluated in [10,11]. Using Clebsch-Gordan (CG) expansion to express Dicke states of N atoms as the product of the states of N − 1 atoms and those of single atom with J ± 1/2 and 1/2, respectively, one identifies for each dissipative mechanism a sum of three terms related to Dicke states with J s = J + s (s = 0, ±1). We simulate the effect of these three terms in the master equation by corresponding quantum jumps of our state vector. Note that the jumps do not correspond to a definite measurement process and outcome but are utilized here merely as a computational tool to simulate the decay of individual atoms in the master equation without embarking into the population of non-symmetric states. In our symmetrical average evolution, the decay of individual atoms thus leads to the transformation of state amplitudes  Table I.
The incoherent excitation (pumping) of individual atoms sketched in Fig. 1 (d) in the main text leads to a similar transformation of state amplitudes implemented with the jump probabilities δp κ,s = δtκ l d J N M |P JM +,s C JM (t) | 2 , where the values of P JM +,s are presented in Table I.
The dephasing of individual atoms retains the value of M and allows changes of J by 0, ±1, as sketched in Fig. 1 (e) in the main text, This process is represented by quantum jumps with the transformation of state amplitudes implemented with the jump probabilities δp d, where the values of P JM z,s are presented in Table I.
To describe the atom loss, we use the CG expansion |JM = j1,j2 m1,m2 JM C j1m1 j2m2 |j 1 m 1 |j 2 m 2 with the CG coefficients JM C j1m1 j2m2 and interpret each term in the expansion as one quantum jump, which removes 2j 1 atoms in the states |j 1 m 1 from the system and leave the 2j 2 = N − 2j 1 atoms in the states |j 2 m 2 . The different jumps will take care of all possible ways of removing the atoms in all possible states. To describe the feeding of new, uncorrelated atoms into the ensemble, we use the inverse CG expansion |j 1 m 1 |JM = j2,m2 j2m2 C j1m1 JM |j 2 m 2 and interpret each term in the expansion as one quantum jump, which adds 2j 1 atoms in the state |j 1 m 1 to the states |JM of N atoms to form the new states |j 2 m 2 of 2j 2 = N + 2j 1 atoms. The different jumps will account for all possible final states by adding atoms in one specified state. Since the probabilities of the jumps are proportional to the squared CG coefficients, the jumps of removing and feeding a single atom are more favored than other jumps. In the case of a single atom loss, i.e. j 1 = 1/2, we will have four different quantum jumps, which remove the single atom in the upper (lower) level and leave the N − 1 atoms in the Dicke states of J ± 1/2 and M − 1/2 (M + 1/2), cf. Fig. 1 (f). These jumps can be implemented with C JtMs (t + δt) = δtγ loss L JM t,s C JM (t) .
Here, we have introduced J t=± = J ± 1/2 and M s=± = M ± 1/2, and the values of L JM t,s are provided in Table I. The jump probabilities are δp l,ts = δtγ loss d J N M |L JM t,s C JM (t) | 2 . In the case of single atom feeding, i.e. j 1 = 1/2, we will also have four different quantum jumps, which add the single atom in the upper (lower) level and leave the N + 1 atoms in the Dicke states of J ± 1/2 and M + 1/2 (M − 1/2), cf. Fig. 1  (g). These jumps can be implemented with where the coefficients F JM t,s are given in Table I. The different processes occur with the jump probabilities δp κ,ts = P s δtκ f eed d J N M |F JM t,s C JM (t) | 2 . For the atom loss we can not control the state of the lost atom and thus should consider all the possibilities. However, for the atom feeding, we can control the jumps by preparing the fed atom in the upper or lower level. For this reason, we introduce the probability P s to describe the fed atom in the upper (s = +) or lower (s = −) level, respectively.

Appendix B: Steady-State Spectrum Calculation
In the main text, we indicate how to calculate the steadystate spectrum with the MCWF method.
In the peculiar case as studied in the main text, the ancillary wave-functions will have only two non-vanishing amplitudes at any point during the stochastic wave-function evaluation. In the evaluation, these amplitudes, which reflect superposition of two adjacent Dicke states, persist despite of the quantum jumps, which ensures that the coherence time can be much longer than the time spent by any typical atom. This guarantees the lasing with millihertz line-width as shown in the main text.