Emergent dark states from superradiant dynamics in multilevel atoms in a cavity

We investigate the collective decay dynamics of atoms with a generic multilevel structure (angular momenta $F\leftrightarrow F'$) coupled to two light modes of different polarization inside a cavity. In contrast to two-level atoms, we find that multilevel atoms can harbour eigenstates that are perfectly dark to cavity decay even within the subspace of permutationally symmetric states (collective Dicke manifold). The dark states arise from destructive interference between different internal transitions and are shown to be entangled. Remarkably, the superradiant decay of multilevel atoms can end up stuck in one of these dark states, where a macroscopic fraction of the atoms remains excited. This opens the door to the preparation of entangled dark states of matter through collective dissipation useful for quantum sensing and quantum simulation. Our predictions should be readily observable in current optical cavity experiments with alkaline-earth atoms or Raman-dressed transitions.


I. INTRODUCTION
Staggering progress in the ability to manipulate quantum matter in atomic, molecular, and optical experiments is opening up new opportunities for quantum sensing, simulation, and computation [1]. So far, most of the effort has been focused on isolating a pair of internal levels (a qubit) to realize paradigmatic models of two-level quantum systems. In recent years, however, experiments have been attaining the degree of control and tunability to study many-body physics using the full atomic multilevel structure. Prominent examples of this direction are experiments working with magnetic atoms [2][3][4][5] featuring strong dipolar interactions and with alkaline-earth(-like) atoms (AEAs), which possess long-lived ground and electronic levels and unique collisional properties [6][7][8][9][10].
The dissipative part of the interactions is responsible for superradiance [42,43], a phenomenon first predicted by Dicke [44], whereby atoms emit light at a collectively enhanced rate due to their collective coupling to a common radiation mode. Superradiant emission is at the heart of laser technologies and has been experimentally observed in a myriad of platforms . In the context of AEAs that feature longer lifetimes than the photons in the cavity, superradiance can give rise to lasers with the highest spectral purity, which could be used to improve the stability of passive atomic clocks [73][74][75][76]. To date, however, the dissipative part of the cavity-mediated interactions has been detrimental for entanglement generation and quantuminformation applications because superradiance strongly reduces the lifetime of the atoms inside the cavity. Experiments have so far circumvented this issue by creating squeezing in the ground-state manifold [13,15,16,25].
Here, by studying the eigenstate structure and cooperative emission properties emerging from the collective coupling of multilevel atoms to a cavity (Fig. 1), we show that superradiance is not necessarily a limitation in multilevel atoms. In contrast to collective (permutationally symmetric) states of two-level atoms, which are all superradiant, we find that multilevel systems support long-lived collective dark states even in the fully symmetric manifold. These dark states are many body in nature and are necessarily entangled. Intuitively, they emerge from the destructive interference of decay channels corresponding to different internal transitions.
These dark states naturally arise as steady states of a superradiant decay process, thus opening the door to the generation of collective entangled states through dissipation [77,78]. This further shows that subradiance can emerge from superradiance in multilevel systems. We note that subradiance has remained experimentally more elusive than superradiance due to its delicate nature, and clear demonstrations of the phenomenon are scarce [14,[79][80][81][82][83][84][85][86][87][88][89][90][91][92]. The dark states studied here are only stable to emission into the relevant cavity modes but are still vulnerable to singleparticle spontaneous emission into free space. Nevertheless, thanks to the distinct separation of timescales over which the dark states are stable, their long-lived character should be experimentally observable in current optical cavity experiments operating with AEAs or with engineered Ramandressed transitions, and even in other settings such as superconducting qubits [93] inside microwave cavities.

A. Background and summary
To put this work into context, it is important to note that, to date, most cavity QED experimental works have focused on engineering effective two-level systems. Only a few setups using alkali atoms have considered cases with three [21,24,27,30,[94][95][96] or more levels [22,[97][98][99], where either the cavity mediates interactions between different ground states or the physics involves at most few atomic excitations. A notable exception is Ref. [71], which studied multilevel aspects of superradiance in the 20-level 1 S 0 -3 P 0 transition of the AEA 87 Sr. Despite this, multilevel cavity physics, especially when atoms have multiple ground and excited states, remains experimentally largely unexplored.
On the theory side, multilevel atoms in cavities have been mostly analyzed in terms of their low-lying excitations [100][101][102], albeit in rather complicated level structures motivated by alkali atoms. More recently, the phases of a multilevel single-mode Dicke model [103] and the collective eigenstates of the vibrational modes of molecules coupled to a cavity [104] have also been investigated. Apart from that, some works have looked at the superradiant transition [105] and squeezing generation in the ground levels [106][107][108]. Thus, except for these scattered works, a surprisingly large knowledge gap exists about the behavior of multilevel atoms inside cavities.
Compared to the few related works described above, the key features and contributions of our work are the following: (i) We consider atoms with degenerate ground (g) and excited (e) levels with associated angular momenta F g ↔ F e [ Fig. 1(c), Sec. II] fully accounting for the corresponding Clebsch-Gordan coefficients. We show that the dynamics is highly sensitive to the underlying multilevel structure.

Atom cavity
FIG. 1. Multilevel atoms in a two-mode cavity. (a) Sketch of the atom-cavity system. The atoms are assumed to be pinned by an external potential (e.g., an optical lattice) and couple to the photon field with strength g c . The photon has two allowed polarizations ⃗ ϵ 1 and ⃗ ϵ 2 , which are orthogonal to each other and to the cavity axis. The photons escape the cavity at a rate κ. (b) After adiabatically eliminating the photons, we obtain an effective atom-only model with photon-mediated collective elastic (χ) and inelastic (Γ) interactions. The inelastic part leads to superradiant decay, which in a multilevel system can generate an entangled steady state. The elastic part is set to zero, χ ¼ 0 for most of the work, except in Secs. IV and XI C. (c) Multilevel internal level structure of the atoms, which is composed of a ground (g) and an excited (e) manifold of states with angular momenta F g and F e , respectively.
(ii) In our system, atoms interact with two cavity modes corresponding to the two different photon polarizations [ Fig. 1(a), Sec. II]. This gives rise to a nontrivial competition between the two decay channels and results in considerably richer physics than two-level atoms. Despite the additional complexity, we find that the problem can be considerably simplified by an appropriate choice of atomic and polarization basis (Sec. II B); this finding is at the heart of this work. (iii) We study cavity-mediated collective decay (Γ) within the collective permutationally symmetric subspace (Sec. III) and generally set cavity-mediated elastic interactions (χ) to zero [ Fig. 1(b)]. This collective regime is naturally realized in cavities, whereas atoms in free space would require short subwavelength interparticle distances for which elastic interactions are important. (iv) We show that, even within the permutationally symmetric manifold, dark states exist and are ubiquitous in l ≥ 6-level systems. These dark states are a result of quantum interference between different internal transitions-a mechanism absent in twolevel atoms. We study the full eigenstate structure and find dark states with up to half of the atoms excited, not just a single excitation (Sec. IV). Moreover, we explicitly show that dark states are entangled (Sec. IV D) and can be populated by superradiance (Sec. V). (v) We introduce a systematic methodology to address multilevel superradiance. One of our main findings is that in many situations the decay is dominated by one polarization, as we justify in Sec. VII. In this case, we find that the system can be described in mean-field approximation in terms of multiple Bloch spheres (Sec. VI A). In this picture, the Bloch-sphere trajectories induced by superradiant dynamics are identical to Rabi drive dynamics (Sec. VI B). The decay rate is given by the (time-dependent) total dipole moment of the atoms, which is a sum over the dipole moments associated with each internal transition. Thus, at the mean-field level, dark states are configurations where the individual dipole moments are nonzero, but the sum cancels out. The mean-field dynamics can be described as the movement of a classical particle on a superradiance potential VðθÞ, which is equivalent to the Rabi excitation profile (Sec. VI C). (vi) We offer a simple and realistic mechanism to prepare (entangled) dark states: Using just a Rabi drive, the atoms can be initialized in a fully coherent initial state, which will naturally evolve through superradiant emission into an entangled dark steady state (Sec. VIII), where a macroscopic fraction of the atoms remains excited. We explicitly show that quantum fluctuations lead to small deviations from the mean-field dark states (minima of the superradiance potential) and are responsible for entanglement generation under superradiance (Secs. V and VIII B). (vii) We also discuss the most general case in which the atoms decay with both polarizations (Sec. IX). In this case, the dynamics is harder to describe, but the concept of superradiance potentials (one for each polarization) can still be used locally, i.e., to determine the stability of a state to decay with either polarization and thus find stable two-polarization dark states. (viii) We propose the use of alkaline-earth(-like) atoms in optical lattices, such as 87 Sr or 171 Yb, as a promising platform to observe the dark states described due to their simple ground-state manifold and the existence of narrow transitions (Sec. X). We show that dark states are surprisingly robust to moderate magnetic fields, inhomogeneous couplings, and the presence of residual elastic interactions χ (Sec. XI). In realistic experimental settings, dark states will inevitably decay due to single-particle spontaneous emission. However, this decay will be slow compared to the collectively enhanced superradiant decay, thus opening a large window of time in which to observe dark states.

II. MULTILEVEL ATOMS IN A CAVITY
A. Two-polarization multilevel model We consider systems composed of multilevel atoms interacting with light inside a cavity (Fig. 1). Atoms are assumed to be pinned by an external confinement potential, and their motion is neglected. There are two main ingredients to this problem: the specific level structure of the atoms, which determines the strength of the allowed transitions, and the properties of the cavity eigenmodes.
We assume that each atom has an internal level structure as shown in Fig. 1(c), which consists of a ground (g) and an excited (e) manifold with total angular momenta F g and F e , respectively. Thus, there are 2F g þ 1 ground and 2F e þ 1 excited states, which can be defined by their angular momentum projection m g and m e onto a fixed quantization axisẑ q . Specifically, we label the ground states as jg m i q , with m ∈ f−F g ; −F g þ 1; …; F g − 1; F g g, and the excited states as je m i q , with m ∈ f−F e ; −F e þ 1; …; F e − 1; F e g. The g and e manifolds are separated by a frequency ω 0 , and all levels within each manifold are assumed to be degenerate.
The atoms are assumed to interact with two degenerate eigenmodes of the cavity, that for now we label asâ 1 andâ 2 . The modes have a frequency ω c detuned from the atomic transition by Δ c ≡ ω 0 − ω c , a linewidth κ, and different orthogonal polarizations ⃗ ϵ 1 and ⃗ ϵ 2 perpendicular to the cavity axis [ Fig. 1(a)]. This is EMERGENT DARK STATES FROM SUPERRADIANT DYNAMICS … PHYS. REV. X 12, 011054 (2022) 011054-3 different from two-level atoms, which can absorb or emit only photons with a single polarization. In the rotating frame of the atoms, the photon dynamics is captured by the HamiltonianĤ c ¼ −Δ c P γ¼1;2â † γâγ and the Lindbladian L κ ðρÞ¼ðκ=2Þ P γ¼1;2 ð2â γρâ † γ −fâ † γâγ ;ρgÞ. The atom-light interaction in dipole approximation (ignoring counterrotating terms) is then given by the multilevel Jaynes-Cummings Hamiltonian which describes the coherent exchange between photons and atomic excitations. Here, 2g c is the single-photon Rabi frequency, which is assumed to be spatially uniform and equal for both polarizations. While the former assumption can be relaxed as we discuss later (see Sec. XI), the latter should be satisfied for an axially symmetric cavity. TheD AE γ operators in Eq. (1) are collective raising or lowering operators that describe the change of the internal state of an atom when it absorbs or emits a cavity photon with polarization ⃗ ϵ γ . Note that, due to the collective coupling, when an atom absorbs a photon, the excitation is coherently shared among all atoms. For two-level atoms, these operators become the usual collective SU(2) spin-raising or -lowering operatorsD AE In the multilevel case, the operatorD AE γ connects all ground states jg n i q and excited states je m i q whose transition dipole moment has a nonvanishing overlap with the photon polarization ⃗ ϵ γ . They are defined whereσ q ab ¼ P i jai q;i hbj q;i is a collective transition operator with i denoting the atom's index. The sum over m, n is assumed to run over all internal levels of e m and g n , respectively.
The vector ⃗ d q e m g n is the spherical part of the transition dipole moment between states jg n i q and je m i q (cf. Ref. [129]) and is given by ⃗ d q e m g n ¼ C m−n n ðê q m−n Þ Ã , where the strength is determined by the Clebsch-Gordan coefficient C p m ≡ hF g ; m; 1; pjF e ; m þ pi and the orientation by the unit vectorsê q m−n , which depend only on the change in angular momentum projection (m − n). These unit vectors are defined with respect to the quantization axis reference frame fx q ;ŷ q ;ẑ q g as π transitions;ê q 0 ≡ẑ q ; ð3Þ Note that only processes with jm − nj ≤ 1 are dipole allowed, i.e., C p m ¼ 0 for p ≥ 2. We work in the limit where the cavity photons do not actively enter in the dynamics, and instead they virtually mediate elastic and dissipative spin interactions between the atoms. In this limit, valid when κ > ffiffiffiffi N p g c or for a fardetuned cavity with jΔ c j > ffiffiffiffi N p g c , we can adiabatically eliminate the photons to obtain an atom-only master LðρÞ ¼ℏΓ X γ¼1;2 where The Hamiltonian part describes photon-mediated collective interactions where two atoms exchange a photon of polarization ⃗ ϵ γ . This term constitutes a multilevel generalization of the collective XX spin model arising in twolevel systems. The two-level model gives rise to a broad class of phenomena such as superconductivity [130], quantum magnetism [19], spin squeezing [131][132][133][134], and dynamical phase transitions [32]. For the multilevel system, what makes the problem particularly complex is that each transition is weighted by a different Clebsch-Gordan coefficient and that, in general, the two operatorsD AE γ for γ ∈ f1; 2g do not commute. The more complex interactions are expected to give rise to rich unitary dynamics awaiting to be explored.
The Lindbladian part describes the collective decay of atoms with either of the two polarizations and is responsible for superradiant emission. While for two-level atoms (D AE γ →Ŝ AE ) the superradiant decay is well known [42], the multiple decay channels and polarizations available in the multilevel case make it hard to anticipate the decay dynamics of the multilevel atoms.
In this work, we study the collective decay dynamics caused by Γ while setting χ ¼ 0. This regime can be easily realized in experiments by setting the cavity on resonance with the atomic transition, i.e., Δ c ¼ 0 [Eq. (7)]. Note that, in this limit, we have Γ ¼ Cγ s , where C ¼ 4g 2 c =ðκγ s Þ is the cooperativity, and γ s is the free-space spontaneous emission rate, which we discuss in Sec. XI D. As we show, the understanding of the dissipative processes can shed light on the rich many-body behavior exhibited by multilevel systems in regimes where analytic insight is possible. We address the role of finite elastic interactions in Sec. XI C.
Before entering into the core discussion of the physics, we discuss next some examples ofD AE γ for specific choices of quantization axis and photon polarizations that we employ throughout this work. This analysis helps provide intuition on what these operators actually do.

B. Multilevel dipole operators: The choice of quantization axis and polarization bases
One of the keys to simplifying multilevel problems is to choose the quantization axisẑ q wisely. Depending on the problem, we employ here one of three different reference frames fx q ;ŷ q ;ẑ q g, q ∈ fV; H; kg, as shown in Fig. 2: the vertical basis V, the horizontal basis H, or the parallel basis k. The latter quantization axis (ẑ k ) is aligned parallel to the cavity axis, whereas the other two (ẑ V ,ẑ H ) are orthogonal to it. Each of these quantization axes defines a basis of atomic states fjg m i V ; je m i V g, fjg m i H ; je m i H g, and fjg m i k ; je m i k g, which can be transformed into a different basis by applying standard rotations; see the Appendix A.
As with the atoms, we are free to choose the polarization basis f⃗ ϵ 1 ; ⃗ ϵ 2 g for the cavity modes fâ 1 ;â 2 g that is best suited to our problem. We employ two different bases as shown in Fig. 2. The first is a decomposition in linearly polarized modes fâ V ;â H g: a mode with vertical polarization ⃗ ϵ V parallel toẑ V , and a mode with horizontal polarization ⃗ ϵ H parallel toẑ H . Alternatively, we can decompose the light modes into their right-and left-handed polarized components fâ R ; Þðâ V þ iâ H Þ. Using these definitions, we define the set of dipole operators associated with the decomposition into linearly polarized modes f⃗ ϵ V ; ⃗ ϵ H g asΠ AE ≡D AE V andΣ AE ≡D AE H .
These operators have different matrix forms depending on the quantization axis used, as depicted in Figs. 2(a) and 2(b). In the vertical basisẑ V , they read [ Fig. 2 In this basis, theΠ AE operator is a sum of π-polarized transitions because the photon polarization ⃗ In contrast, theΣ AE operator is a sum of σ AE -polarized transitions, because the photon polarization ⃗ ϵ H ¼ŷ V has a 1= ffiffi ffi 2 p overlap with the dipole of both Thus, if an atom is in the state jg m i V and absorbs an ⃗ ϵ V -polarized photon, it will be excited to je m i V . However, if the atom absorbs an ⃗ ϵ H -polarized photon, it will be excited to a superposition of je mAE1 i V . In the horizontal basisẑ H , the role of these operators is reversed as shown in Fig. 2(b); i.e., Π AE is a sum of σ AE -polarized transitions, whereasΣ AE is a sum of π-polarized transitions.
Analogously, we define the set of operators associated with the decomposition into circularly polarized modes f⃗ ϵ L ; ⃗ ϵ R g asL AE ≡D AE L andR AE ≡D AE R . These operators acquire their most simple representation when expressed in the quantization axis parallel to the cavity axisẑ k aŝ FIG. 2. Atomic bases and quantization axes. The panels exhibit some properties for different choices of atomic basis: (a) q ¼ V, Each panel shows (i) the frame of reference fx q ;ŷ q ;ẑ q g with quantization axisẑ q , (ii) the typical cavityphoton-polarization decomposition that we use for each choice of atomic basis, and (iii) a pictorial representation of the dipole operators associated with each of these polarizations when written in the q atomic basis. Specifically, each dipole operator is a (weighted) sum over allσ q αβ operators for which the pair of states jαi q and jβi q is connected by a colored arrow. The operators in panel (a) correspond to Eqs. (8) and (9) and the operators in panel (c) to Eqs. (10) and (11).

III. DYNAMICS IN PERMUTATIONALLY SYMMETRIC MANIFOLD
In this paper, we are interested in studying the collective decay dynamics of the model of Eqs. (5) and (6). Since the master equation is written in terms of collective operators only, the dynamics takes place within the Hilbert subspace of permutationally symmetric (PS) states. PS states are those which are invariant to the exchange of the internal state of any two atoms; i.e., all atoms are essentially indistinguishable, and any information is equally spread among all of them.
Each of our multilevel atoms has l ¼ 2ðF g þ F e Þ þ 2 internal levels. For l-level atoms, a convenient basis of PS states can be defined by l-tuples ⃗ n ¼ ðn 1 ; …; n l Þ of nonnegative integers specifying how many atoms are in each of the l levels. In order to reflect the ground and excited internal level structure of the atoms, we denote these PS basis states as j⃗ ni q ≡ n e −F e ; …; n e F e n g −F g ; …; n g F g q ; ð12Þ where n α ∈ N gives the number of atoms in state jαi q , and P m n g m þ P m n e m ¼ N. The state of Eq. (12) is essentially the normalized sum over all ð N ⃗ n Þ ¼ N!=ð Q α n α !Þ distinct permutations of states with N atoms in the levels specified by the n α integers. Because of the permutational symmetry, collective operators act asσ q αα j⃗ ni q ¼ n α j⃗ ni q and where α ≠ β andα ¼ ð0; …; 1; …; 0Þ is a unit vector with a 1 at the αth position. The size of the permutationally symmetric Hilbert space H l sym is given by the number of distinct sets of l integers n α that add up to N, which is given by the binomial coefficient Equation (14) highlights the inherent complexity of multilevel systems. Even in the PS subspace, the collective Hilbert space of atoms with l internal levels scales as N l−1 .
Only for two-level atoms the scaling is linear with N. Thus, although the collective nature of the problem considerably simplifies its description, the multiple internal levels still constitute a challenge. While we can use exact diagonalization (ED) methods for small N on the order of tens (size ofρ scales as N 2ðl−1Þ ), we employ different approximations to study large N: a mean-field (MF) approximation, a truncated Wigner approximation (TWA), and a cumulant expansion to second order (cumulant). These methods are outlined in Appendix B.

IV. SYMMETRIC EIGENSTATES
In order to understand the decay dynamics of our multilevel model of Eqs. (5) and (6), it is extremely useful to look at the system's eigenstates within the collective manifold. Specifically, this provides information about the lifetime of many-body eigenstates and the states they decay into.
The problem of finding the multilevel eigenstates considerably simplifies when we employ the basis defined by the parallel quantization axisẑ k together with the decomposition into right and left operators,L AE andR AE . Thus, we consider the effective non-Hermitian Hamiltonian [obtained as usual from Eqs. (5) and (6) without the recycling term] given bŷ As we mentioned above, we set χ ¼ 0 in the following sections by choosing to work on resonance, Δ c ¼ 0 [Eq. (7)]. However, note that the form of the eigenstates ofĤ eff is independent of the specific values of χ and Γ. In particular, eigenstates that are dark to decay (Sec. IV C), which are at the core of this paper, remain dark for χ ≠ 0. A nonzero χ does, however, modify the dynamics, as we discuss in Sec. XI C. The Hamiltonian (15) conserves the total number of excitationsN e ¼ P mσ q e m e m , so its eigenstates can be labeled by the number of atoms N e present in the excited states. The (right) eigenstates jψ k i ofĤ eff have complex eigenvalues,Ĥ eff jψ Here, ℏε k ∝ ℏχ is the energy shift of the eigenstate, and γ k ∝ Γ represents its decay rate into states with one excitation less. Note that energy shifts and decay rates are simply related by ε k =χ ¼ γ k =Γ.
The effective Hamiltonian Eq. (15) has a number of conservation laws which become obvious when usingL AE , R AE , and theẑ k basis. This allows us to reduce the size of and sometimes analytically solve the eigenvalue problem. In the following subsections, we discuss specific two-, four-, and six-level examples in order to illustrate some of these conservation laws, the rich variety of eigenstates available, and how they determine the decay dynamics of multilevel atoms. For other cases, we refer to Appendix C.

A. Two-level Dicke states
In the two-level limit, whereĤ eff ¼ ℏðχ − iΓ=2ÞŜ þŜ− , there is only one possible polarization and decay channel. In this case, the N þ 1 collective eigenstates are the wellknown Dicke states [44,135], which simply correspond to the PS states j k N−k i, where k excitations are symmetrically shared among all atoms. What is important to note is that there is only one Dicke state for each fixed value of k. Because of this, there is only one possible decay path: The atoms collectively decay down the ladder of Dicke states (k → k − 1 → k − 2 → …) all the way down to the ground state k ¼ 0 [see Fig. 3(a)]. The ensuing dynamics is determined by the decay rate of the Dicke states, which is given by The ðN − k þ 1Þ prefactor is responsible for the enhancement with respect to the decay of independent atoms (γ k =Γ ¼ k), which is the hallmark of superradiance. As we will see in a moment, the fact that all two-level collective states are superradiant does not generalize to other level structures.

B. Four-level eigenstates: Dicke generalization
Four-level systems [l ¼ 2ðF g þ F e Þ þ 2 ¼ 4] are slightly more complicated than two-level systems because they can decay via two polarizations, i.e., viaL − orR − . Nevertheless, in four-level systems the operatorsL − andR − have the simplifying property that each of them acts only on a single atomic transition in theẑ k atomic basis. Because of this, the four-level eigenstates ofĤ eff , Eq. (15), are simply the PS states of Eq. (12). These eigenstates can be seen as a multilevel generalization of Dicke states.
To explore the properties of four-level systems, we focus here on the case with ðF g ; F e Þ ¼ ð1=2; 1=2Þ. In this level structure,L AE connects the states fjg 1=2 i k ; je −1=2 i k g, andR AE connects the other two states, fjg −1=2 i k ; je 1=2 i k g. Thus, we can parametrize the multilevel Dicke states as where The numbers N R and N L are conserved quantities and correspond to the number of atoms in the levels addressed byR AE andL AE , respectively (cf. Appendix C 1). These four-level Dicke states, Eq. (17), decay with a rate where for simplicity we write γ k L ;N L ;k R ;N R → γ. Up to the 2=3 prefactor coming from the Clebsch-Gordan coefficient, this expression is essentially the sum over two two-level decay rates corresponding to theL − andR − decay operators; cf. Eq. (16). This is becauseL AE andR AE commute for this level structure. Despite this, the four-level system cannot be simply viewed as two independent two-level systems.
In particular, the spectrum of decay rates is richer than for two-level atoms, as exemplified for N ¼ 10 and different values of k ¼ k L þ k R in Figs. 3(b) and 3(c). The distribution of decay rates spans from single-particle decay (γ=Γ ∼ k) all the way up to N-enhanced superradiant decay (γ=Γ ∼ Nk). Broadly speaking, the more atoms are present in a state where an excitation can decay, the larger the superradiant enhancement. The states with the minimal decay rate are of the form j k N−k 0 0 i k , whereas the maximally decaying states look like j k 0 0 N−k i k . This shows two key features that we will encounter repeatedly: (i) Multilevel systems can harbor highly entangled states with small We show results for a two-level system, a four-level system with F g ¼ F e ¼ 1=2, and a six-level system with F g ¼ 1=2 and F e ¼ 3=2.
decay rates, and (ii) the most superradiant states are two level in nature. As opposed to two-level atoms, the decay path of four-level atoms is not unique. Every time an atom decays, there is a bifurcation [see Fig. 3(a)]: A state with excitations ðk L ; k R Þ can decay to ðk L − 1; k R Þ or ðk L ; k R − 1Þ viaL − orR − , respectively. This leads to a complicated cascade of decay paths weighted with different probabilities. However, note that sinceL AE andR AE commute, for an initial state such as Eq. (17) with fixed N R and N L all decay paths eventually converge to the same ground state, j 0 N R 0 N L i k . This allows us to gain analytical insight into the properties of the final distribution of ground states resulting from the decay, which we further explore in Appendix G.

C. Six-level eigenstates: Beyond Dicke and dark states
Multilevel systems with l ≥ 6 are more complex due to the possibility of interference between different internal levels. This can be seen from the fact that the operatorsL AE andR AE now act on multiple atomic transitions. As a consequence, the PS states of Eq. (12) are no longer eigenstates of the system. Instead, eigenstates are generally given by superpositions of multilevel Dicke states with coefficients that depend on Clebsch-Gordan coefficients. Interestingly, this allows for the appearance of dark collective eigenstates, which we define as states with a nonzero amount of excitations yet zero decay rate due to destructive interference of decay channels. Recall that "dark" refers to decay into the cavity modes; free space emission is discussed in Sec. XI.
To illustrate these features, we consider the six-level system ðF g ; F e Þ ¼ ð1=2; 3=2Þ. Figures 3(b) and 3(c) show the distribution of decay rates for this six-level system with N ¼ 10 and different numbers of excitations k. Remarkably, we find dark states with up to k ∼ N=2 excitations. The remaining nondark eigenstates again span a wide range of decay rates with the maximal decay rate corresponding to the two-level limit of Eq. (16). These maximally superradiant states correspond to states such as j 0 0 0 0 N−k k i k , which again involve only two levels. To gain some insight into the mechanism behind dark states, it is useful to consider a specific example. For k ¼ 1 excitations, we find that the following (unnormalized) states are dark eigenstates: where N A ; N B ≥ 1, and give the number of atoms in the fg −1=2 ; e −3=2 ; e 1=2 g (fg 1=2 ; e −1=2 ; e 3=2 g) levels, and they are dynamically conserved; see Appendix C 1. It is straightforward to show using Eq. (13) that the states of Eq. (19) are zero eigenstates of both decay operators, which are given The darkness of the states in Eq. (19) sensitively depends on a careful cancellation of the decay processes e 3=2 → g 1=2 and e 1=2 → g −1=2 made possible by the choice of amplitudes. To illustrate this, note that the (unnormalized) state is also an eigenstate ofĤ eff but with a superradiant decay rate γ=Γ ¼ N B þ N A =3 due to constructive interference. Finding analytical formulas for the eigenstates with k ≥ 2 becomes increasingly hard. Nevertheless, we can find expressions for dark states with arbitrary k in the case where the atoms can decay only with, e.g.,R − , i.e., when the population of the je −3=2 i k and je −1=2 i k states is zero. Under this restriction, all dark eigenstates that exist are given by The relationship between the coefficients a N A r follows from the fact that the decay channels of the different PS states in Eq. (21) are canceled in pairs. There is exactly one eigenstate for each allowed k and N A , which gives a total of P bN=2c k¼1 ðN − 2k þ 1Þ dark states. Analogous expressions apply to states that can decay only withL − .
Finding general expressions for the remaining dark and bright eigenstates is involved. However, the eigenvalue problem can be considerably simplified using conservation laws as detailed in Appendix C. In general, we find that dark states are ubiquitous in multilevel atoms with l ≥ 6 internal levels.

D. Dark states are entangled
An important property of dark states is that, as long as the atomic structure does not admit trivial single-particle dark states [which is only the case for the ðF g ; F e Þ ¼ ð0; 1Þ level structure], they are necessarily entangled. To prove this, one can show (see Appendix D) that product states of N atoms cannot be dark because they necessarily fulfill for at least one of the polarizations γ. This implies that many-body dark states need to be entangled. Because of this property, we sometimes refer to dark states of the full quantum theory as quantum dark states.
To illustrate the entanglement content of quantum dark states, we show in Fig. 4 the Renyi entropy R 1 ≡ − log 4 ρ 2 1 , where ρ 1 is the single-particle reduced density matrix for each of the dark eigenstates of Eq. (21) with different k and N A . All dark eigenstates shown are entangled since R 1 > 0. Two-level entangled states fulfill 0 < R 1 ≤ 0.5, which is saturated by the Dicke state j N=2 N=2 i. However, the dark states of Eq. (21) are effectively four level and fulfill instead 0 < R 1 ≤ 1. While none of the dark eigenstates shown saturates the bound, some of them are rather close to R 1 ¼ 1.

V. DARK-STATE ENTANGLEMENT FROM COLLECTIVE DECAY
The existence of dark states raises a question: What is the role of such dark states in the decay dynamics of multilevel atoms? At the level of eigenstates, we can already anticipate that atoms might decay into and get stuck in a dark state [ Fig. 3(a)]. For example, for N ¼ 2 one can easily show that the doubly excited state j 0 0 0 1 0 1 i k decays partially into the singly excited dark state of Eq. (19). However, it is a priori not clear how this behavior depends on the initial state and N.
As a teaser, we show in Fig. 5 the ED decay dynamics of the total excited-state population n e ¼ ð1=NÞ P m hσ q e m e m i for a six-level system with F g ¼ 1=2 and F e ¼ 3=2. The system starts in a product state with each atom in a superposition ðjg −1=2 i k − jg 1=2 i k Þ= ffiffi ffi 2 p , it is then excited by a laser of pulse area θ 0 and polarization ⃗ ϵ R (see caption and Sec. VI for details), and is finally let to decay with Eq. (6). In Figs. 5(a) and 5(b), we plot the dynamics for the cases θ 0 ¼ 1π and 1.3π, which shows that for N > 1 the system ends up in a (quantum) dark state with nonvanishing excitations as t → ∞. Following the discussion of the previous section [see Eq. (22)], this (quantum) dark steady state of the collective decay dynamics is necessarily entangled.
This behavior can be qualitatively understood from the eigenstate structure discussed in Sec. IV C. Right after excitation, the system is in an unentangled superposition of different j 0 0 Each of these states can be written as a superposition of bright eigenstates and the dark eigenstate of Eq. (21). For each k and N A the population in the dark eigenstate remains stuck, whereas the various bright eigenstates will decay to a superposition of bright and dark eigenstates with k − 1 excitations and same N A . This process repeats itself until where ρ 1 denotes the single-particle reduced density matrix computed for dark states of a six-level system with F g ¼ 1=2 and F e ¼ 3=2. We show all the dark states for N ¼ 50 with excitation number k ∈ ½1; N=2 and number parameter N A which have no occupation in the je −3=2 i k and je −1=2 i k states as defined in Eq. (21). All dark states shown fulfill R 1 > 0 and are, therefore, entangled.
and is excited with a expð−iτĤ Ω Þ pulse, whereĤ Ω ¼ ðΩ=2ÞðR þ þR − Þ and the pulse area is given by θ 0 ≡ jΩjτ. the final steady state is reached, which is an entangled mixed state of multiple dark and ground states.
The properties of the final state sensitively depend on the initial rotation θ 0 . Figures 5(a) and 5(b) show that as N is increased, the final fraction of excitations seems to approach zero for θ 0 ¼ 1π but a nonzero value for θ 0 ¼ 1.3π. In Fig. 5(c), we confirm this behavior by plotting the long-term value of n e as a function of θ 0 . Intriguingly, we find signatures of a potential transition between different values of n e as N increases.
The large-N behavior of Fig. 5 cannot be explored using ED due to memory constraints. Moreover, even if simulations were feasible, quantitatively understanding the dynamics in terms of the eigenstates would be complicated due to the complexity of the eigenstate structure. Therefore, we study in the next section the large-N limit using a meanfield approximation. The physical scenario of Fig. 5 is revisited in Sec. VIII.

A. Single-mode model with multiple Bloch spheres
To shed light on the decay dynamics of multilevel atoms for large N, we first consider the case where the decay is dominated by only one of the two polarization modes. In such cases, one can describe the system in terms of a singlemode model where the multilevel atoms are coupled to a single cavity mode of polarization ⃗ ϵ, leading to [cf. Eq. (6)] Here, the polarization ⃗ ϵ and the operatorsD AE are related through Eq. (2) and could take any of the forms introduced above.
The key to understanding this model is the realization that we can always find an atomic basis fjgðαÞi; α ∈ ½1; 2F g þ 1g and fjeðαÞi; α ∈ ½1; 2F e þ 1g such that the operatorD AE couples the ground and excited levels in distinct pairs feðαÞ; gðαÞg with no overlap between pairs. We write it generically asD − ¼ ðD þ Þ † , where the c α ∈ R are related to Clebsch-Gordan coefficients, and P α is assumed to run from 1 to 2 minðF g ; F e Þ þ 1. Note that some of the internal levels may not be coupled byD AE [e.g., in Eqs. (10) and (11)]. This expression (24) is essentially a sum of two-level raising or lowering (Pauli) operators acting on disjoint twolevel subspaces of the multilevel atom. Hence, we refer to it as the multi-two-level form ofD AE .
A general way to write anyD AE in the form of Eq. (24) is as follows. If the mode considered has a linear polarization ⃗ ϵ, then we choose an atomic basis with the quantization axiŝ z q parallel to ⃗ ϵ. As shown in Fig. 2,Π AE takes this form when written in theẑ V atomic basis [see Eq. (8)], and Σ AE in theẑ H basis. For circular polarization ⃗ ϵ AE , the operatorsL AE andR AE have this form in theẑ k basis [see Eqs. (10) and (11)]. For the more generic case of elliptical polarization, we refer to Appendix F.
We analyze the superradiant decay dynamics of this single-mode model using a MF approximation. The beauty of the transformation into multi-two-level form, Eq. (24), is that it allows us to describe the system in terms of multiple Bloch spheres, each one associated with a different twolevel subspace fjgðαÞi; jeðαÞig; see Fig. 6(a). Specifically, we can define a set of α-spin variables for each of the above pairs of levels as withσ þ α ¼σ eðαÞgðαÞ andσ z α ¼σ eðαÞeðαÞ −σ gðαÞgðαÞ . As we see, the length of each of these α-spins is conserved by the mean-field dynamics, Hence, we can describe these variables by vectors ⃗ S α ¼ ðS x α ; S y α ; S z α Þ on the surface of α-Bloch spheres of fixed radius s α . While in a two-level system we have s α ¼ N=2, in a multilevel system s α is given by half the total number of atoms in the fgðαÞ; eðαÞg subspace. Thus, s α can be anything between 0 and N=2, with the total sum being P α s α ¼ N=2. In principle, there are additional single-body observables which involve coherences between the α-Bloch spheres. However, these observables decouple in MF from the ⃗ S α variables and will be ignored henceforth.

B. Mean-field decay dynamics
We consider an excitation-decay scenario as depicted in Fig. 6(b). The atoms are all initially prepared in some arbitrary ground state jgi. This means that all α-Bloch vectors start in the south pole S z α ¼ −s α . The atoms are then coherently excited through the cavity using a laser pulse of duration τ and on resonance with the atomic transition, which is described by the Rabi Hamiltonian where Ω ¼ Ω x − iΩ y . This leads to a rotation of each αspin as described by the MF equation of motion where the torque vector ⃗ Ω ¼ ðΩ x ; Ω y ; 0Þ is the axis of rotation. After preparation, the atoms are left to decay freely through Eq. (23). This is described by the MF equation of motion Here, the torque vector ⃗ D ⊥ ¼ ð−D y ; D x ; 0Þ is the vector perpendicular to the total dipole moment ⃗ D ¼ ðD x ; D y ; 0Þ defined by D x;y ¼ P α c α S x;y α . Note that, as anticipated, the total spin length s α of each of the α-Bloch vectors is conserved by both Eqs. (28) and (29), and that the value of s α is determined by the choice of initial ground state jgi: Every state with vanishing dipole moment ⃗ D constitutes a stationary fixed point of the MF decay dynamics, Eq. (29). In a two-level system, the torque vector becomes ⃗ D ⊥ → ð−S y ; S x ; 0Þ, which can vanish only at the north and south poles of the Bloch sphere. However, in a multilevel system all α-spin vectors are coupled to each other through a common torque field ⃗ D ⊥ which is given by the sum over the dipole moments of all the α-spins. Because of this, we can have a vanishing total dipole moment where the individual dipole moments of the α-spins are nonzero but cancel each other out. This is the main ingredient for the existence of multilevel dark states at the MF level.
Two observations are central in our analysis of the decay dynamics. First, the equations of motion for Rabi oscillations and for collective decay, Eqs. (28) and (29), are identical if one identifies ⃗ Ω ↔ Γ ⃗ D ⊥ . Second, the direction of the torque vector for the decay dynamics is conserved, Since ⃗ Ω is trivially constant too, this implies that both Rabi excitation and collective decay follow the same trajectories on the surface of the α-Bloch spheres.
Equation (30) implies that the dynamics of each ⃗ S α happens on a fixed plane perpendicular to ⃗ D ⊥ . Because of this, we can define spin components S k α and S ⊥ α which are parallel and perpendicular to ⃗ D ⊥ , respectively, such that S k α ¼ const. Thus, the dynamics of the remaining variables S ⊥ α and S z α happens on a circle of radius [cf. Eq. (26)] Since the MF equations for different α, Eqs. (28) and (29), differ only in the α-dependent rotation speed c α , we conclude that we can parametrize all Bloch spheres in terms of one parameter θðtÞ as [ Fig. 6 Inserting this into Eq. (29) leads to Thus, the single-mode superradiance problem for multilevel atoms has been reduced to a single equation for the angle variable θðtÞ. The values of r α and φ α are determined by the state at t ¼ 0. In the excitation-decay scenario described above, we have that Further details on this derivation are given in Appendix E.

C. Superradiance potential
The dynamics of θðtÞ given by Eq. (33) can be conveniently described in terms of a superradiance potential VðθÞ as (independent from our work, we found that the possibility of defining a superradiance potential is briefly mentioned for a very specific case study in the last chapter of Ref. [42]) where Note that VðθÞ implicitly depends on the initial state through the parameters r α and φ α (we suppress this dependence for simplicity).
The dynamics of θðtÞ can then be understood as the movement of a classical particle in the potential VðθÞ in the limit of infinite friction; see Fig. 6(c). The dynamics starts at θð0Þ ¼ θ 0 , and the slope of the potential determines in which direction the classical particle θ will move. When the slope of VðθÞ is positive the particle moves to the left, and when it is negative the particle moves to the right.
The extrema of VðθÞ represent stationary states of the MF dynamics, which are states where the average dipole moment of the atoms is zero, D ¼ 0. We call these states MF dark whenever their excited-state population is nonzero in order to distinguish them from the quantum dark states discussed in Sec. IV D. Since MF dark states are unentangled, they cannot be perfectly dark in the full quantum theory due to Eq. (22). The relationship between quantum and MF dark states is explored further below.
The stability of these MF stationary points to (quantum) fluctuations depends on the curvature of the potential around them, ðd 2 =dθ 2 ÞVðθÞ. A stationary point θ SS is stable if the curvature is positive (maximum) and unstable if negative (minimum). If the curvature is zero, higher orders need to be computed. Saddle points represent a curious case because at the MF level they are stable to perturbations in one direction but unstable in the other direction.
Because of the similarity of Eqs. (28) and (29), the superradiance potential is directly related to the dynamics of the total excited-state population (or inversion) under pure Rabi oscillations through Here, S z α ðτÞ is obtained by starting from a state specified by the same r α and φ α as in Eq. (36) and then evolving with Eq. (28) (see Appendix E). The above expression is useful for computing VðθÞ. According to it, the states corresponding to the minima of the Rabi oscillation profile are stable MF dark states, maxima are unstable MF dark states, and saddle points require further investigation. In a two-level system, the potential is VðθÞ − 1=2 ∼ cos θ, and the minima and maxima correspond to the (stable) south and (unstable) north poles of the Bloch sphere.

VII. CLEBSCH-GORDAN COEFFICIENTS VS COHERENCES
Before delving into the phenomena predicted by the single-polarization model of the previous section, it is useful to understand better why this approximation might be valid. As we discuss in Sec. IV, multilevel atoms can decay via multiple paths. The probability of taking one path or another depends on the choice of initial state as well as Clebsch-Gordan coefficients, which we explore in the following.
Following Sec. VI, the superradiant decay of two-level systems [42] can be described by a single Bloch vector ⃗ S ¼ ðS x ; S y ; S z Þ of radius s ¼ N=2 and the superradiance potential VðθÞ ¼ 1 2 ð1 − cos θÞ. If we start close to the north pole θ 0 ≈ π and follow the potential, the collective dipole moment of the atoms [see Eq. (29)], which is proportional to the atomic coherences jS ⊥ j, first increases until it becomes maximal (approximately N=2) at the equator, and after that it decreases until it reaches zero at the south pole. The same behavior is displayed by the emitted light intensity, which follows one to one the decay rate of the atoms-both are proportional to ½S ⊥ ðtÞ 2 . This leads to the emission of a coherent superradiant pulse of light whose peak intensity scales as approximately N 2 =4 at the equator.
Since initially at the north pole there are no coherences, it takes the system time to build them up. This can be quantified by the delay time t D that the Bloch vector needs to reach the equator. For θ 0 ≠ π, the MF equations can be solved exactly [42] and lead to S z ðtÞ ¼ where s 0 z ≡ S z ð0Þ=N. The closer to the north pole, the longer it takes for superradiance to build up. For θ ¼ π, the MF approximation predicts an infinite delay time. However, this is a breakdown of the MF approximation, since at the north pole the dynamics is driven by quantum fluctuations. An estimate of t D can be obtained by noticing that for θ ∼ π, For a coherent state, we have hjŝ þ 0 j 2 i ¼ 1=N, which is consistent with a delay time given by [42] which grows logarithmically as a function of N (see also Sec. IX E). Similar considerations dictate the behavior of superradiance in multilevel systems. As a simple multilevel toy model, we consider first superradiance in an effective Λ-type three-level system with one excited level and two ground states [see Fig. 7(a)]. Such a system can be realized using strong magnetic fields to make transitions to other levels off resonant. As an example, we consider the jg AE1=2 i V and je −1=2 i V states of the ðF g ; F e Þ ¼ ð1=2; 3=2Þ level structure, which we denote as g AE1=2 ≡ g AE and e −1=2 ≡ e for simplicity. Note that we work in this section in the verticalẑ V atomic basis. The e state can decay viaΠ − into g − or viaΣ − into g þ , with single-atom decay probabilities Thus, theΠ − decay into g − is 4 times faster than theΣ − into g þ . If we prepare an initial state with all atoms in cosðθ 0 =2Þjg þ i V þ sinðθ 0 =2Þjei V , there are two competing behaviors one needs to consider. On the one hand, the coherences between g þ and e will trigger superradiance in that channel; on the other hand, Clebsch-Gordan coefficients make decay into g − the preferred option. Moreover, every time an atom decays into one ground state, the probability that the next atom decays into the same ground state increases (see also Ref. [136] and Appendix G).
An investigation of this competition is shown in Fig. 7(b), which shows the final (normalized) population difference between the two ground states n g þ − n g − as a function of the tilt angle θ 0 computed with TWA for different N. When the initial coherences are large (θ 0 → π=2), coherences win and most atoms decay into g þ ; however, for small initial coherences (θ 0 → π), the stronger transition to g − dominates. Remarkably though, if we keep θ 0 fixed and increase N, the initial coherences eventually always win and most atoms end up in g þ . This is because in the absence of coherences, superradiance is delayed by a factor log N [Eq. (39)].
This analysis illustrates the fact that multilevel systems can predominantly decay through one channel even in the presence of two decay channels. Therefore, it justifies the relevance of the single-mode model studied in Sec. VI.

VIII. SINGLE-POLARIZATION DARK STATES
In this section, we use numerical simulations to benchmark the simple mechanism provided by the superradiance potential VðθÞ for the emergence of not only MF dark states but also quantum dark states. Surprisingly, we find that in regimes where VðθÞ predicts the formation of a MF dark state, the 1=N beyond MF corrections induced by quantum  FIG. 7. Clebsch-Gordan coefficients vs initial coherences. (a) Three-level toy model with two ground states g AE and one excited state e. The single-particle decay rates p e→g − ¼ 0.8 and p e→g þ ¼ 0.2 correspond to the jg AE1=2 i V and je −1=2 i V states of the ðF g ; F e Þ ¼ ð1=2; 3=2Þ level structure. The initial state is a superposition of jei and jg þ i, which we depict through the corresponding ðe; g − Þ and ðe; g þ Þ (normalized) Bloch vectors. Only the latter Bloch sphere has nonzero initial coherences captured by the polar angle θ 0 . (b) Population imbalance n g þ − n g − between the two ground states obtained as t → ∞ as a function of the initial tilt θ 0 for different atom numbers N computed with TWA for 10 4 trajectories.
fluctuations generate small but necessary correlations as the system approaches a quantum dark state, which happens to be rather close to the MF dark state. We test in this way the validity of the MF picture to capture the underlying physics of the quantum system.
In the following, we consider scenarios where the atoms start in some ground state jg m i q defined with respect to a quantization vectorẑ q but are Rabi excited using a laser drive with a polarization ⃗ ϵ that is not parallel toẑ q . Thus, from the point of view of the drive, the initial ground state will look like a superposition of different Zeeman states. Alternatively, such superpositions can also be achieved using magnetic fields or microwave drives to perform rotations within the ground-state manifold.

A. Mean-field numerical results
We consider here the decay dynamics of six-level atoms with F g ¼ 1=2 and F e ¼ 3=2 (Fig. 8). We choose this level structure because it allows us to work in a regime where only one polarization is relevant and because nontrivial collective dark states exist only for atoms with l ≥ 6 internal levels; see Sec. IV. We assume that the atoms are initially prepared in the state jg −1=2 i V defined with respect to the verticalẑ V basis and are then Rabi excited with a laser of right-handed polarization ⃗ ϵ R and pulse area θ 0 ¼ jΩjτ (same configuration as in Sec. V and Fig. 5). The associated right operator has the multi-twolevel form of Eq. (24) when written in the parallelẑ k basis and is given byR Importantly, theR AE operator couples only to the je m>0 i k excited states. Because of this, the left operatorL AE , which couples only to the je m<0 i k states, completely decouples from the dynamics. Thus, this becomes exactly a singlepolarization problem taking place in an effective four-level subspace. We get one Bloch sphere for each of the pairs of states fjg m i k ; je mþ1 i k g with m ¼ AE1=2. The radius of each Bloch sphere is determined by the initial population of each of the two-level subspaces. Writing the initial ground state in thê z k basis leads to which implies that both Bloch spheres have radius N=4.
In the MF approximation, the atoms start at the position θ ¼ θ 0 of the potential and start rolling down as dictated by the slope. At long times t → ∞, the atoms settle down in the first minimum that they encounter. The superradiance potential for this scenario is given by and is plotted in Fig. 8(a) in green. This potential shows a plethora of minima, i.e., stable MF dark states. Because of the ffiffi ffi 3 p incommensurable ratio between the two Clebsch-Gordan coefficients in Eq. (41), the Rabi oscillations never fully rephase, and the number of distinct MF dark states in this simple situation is infinite. In Fig. 8(b), we show the (normalized) excited-state population n e ðt → ∞Þ that remains at long times as a function of the initial pulse θ 0 . The system shows sharp transitions between different long-time MF dark states wherever VðθÞ has a maximum. This is because at the MF level an initial infinitesimal displacement to the left or right of a maximum is sufficient to make the system decay into the minimum on the left or right of the maximum. show the time evolution of each Bloch vector for two slightly different initial pulse areas θ 0 ¼ 1π; 1.3π. For θ 0 ¼ 1π, the atoms simply decay back to the south pole, which corresponds to the ground state jg −1=2 i V initially prepared. However, for θ 0 ¼ 1.3π the atoms end up in a MF dark state. The Bloch vectors rotate in opposite directions in these two cases, reflecting the different sign of the potential slope. Interestingly, even though the total excited-state population continuously decreases, for θ 0 ¼ 1.3π the individual population of je 1=2 i k or je 3=2 i k can in principle grow. This type of incoherent population exchange between the excited levels is due to the anticommutator terms in Eq. (6) and can be observed in two-level systems outside the collective manifold as well.

B. Quantum fluctuations
We now explore what happens to the MF dark states when quantum fluctuations and beyond-mean-field effects are included. In Fig. 8(b), we plot the t → ∞ excited-state population for different values of N obtained using TWA (color points). As N → ∞ the TWA results perfectly converge toward the MF prediction (black line). In fact, the full time evolution of the populations is perfectly captured by MF in this N → ∞ limit, as long as the initial state is not too close to a transition point.
At a transition point, i.e., when the atoms start at a maximum of VðθÞ, the mean dipole moment is zero, and hence, the decay is solely driven by quantum fluctuations. This is analogous to the case of a two-level system starting at the north pole of the Bloch sphere. The main difference is that, in our multilevel case, quantum fluctuations can make the atoms decay into different dark states depending on the direction of the fluctuation. For finite N, these quantum fluctuations lead to a smoothening of the transitions as shown in the TWA results of the average excitation fraction in Fig. 8(b). (We note that TWA can lead to unphysical negative populations which, however, converge to zero as N → ∞. Moreover, we checked that cumulant results also converge toward MF as N → ∞, although the results deviate from TWA close to the transition point for finite N.) This is again consistent with the ED results of Fig. 5(c). To shed more light on the case where the atoms start at a transition point, we study the full distribution function of the excitations in Appendix H.
These results indicate that stable MF dark states become asymptotically quantum dark states as N → ∞. However, for any finite N, MF dark states cannot be identical to quantum dark states because the latter must be entangled, as we discuss in Sec. IV D. Thus, if we start at a minimum of VðθÞ, i.e., at a MF dark state, atoms will necessarily decay and emit light. More specifically, atoms starting at a MF dark state will have ffiffiffiffi N p fluctuations leading to dipole moment fluctuations of order hR þR− i ∼ N; cf. Eq. (22). For a state to be truly dark, these fluctuations need to decay. Therefore, the decay of a MF dark state into a quantum dark state happens through the emission of photons with order N light intensity, which leaves the atoms slightly entangled. Note that this contrasts with typical superradiant emission which is of order N 2 .
To illustrate this behavior, we show in Fig. 9 the evolution of the system computed with cumulant starting at the θ 0 ¼ 2.123π minimum; cf. Fig. 8(a). Figure 9(a) shows the decay of the total excited-state population with respect to its initial value, n e ðtÞ − n e ð0Þ. The lines for different N [see legend in Fig. 9(b)] reveal a population decay of order 1=N. Figure 9(b) shows the intensity of the light emission into the right-handed polarized mode normalized as hR þR− i=N. The collapse of the results for different N confirms the above picture. Interestingly, this shows that a MF dark state can become a quantum dark state through the decay of just Oð1Þ number of photons.

A. Orthogonal superradiance potential
In order to shed light on the two-polarization problem, it is useful to extend the single-mode toolkit developed in Sec. VI. We consider a situation where the atoms are excited (through the cavity) via a Rabi driveD AE 1 of some polarization ⃗ ϵ 1 , but the decay happens with bothD − 1 as well asD − 2 of orthogonal polarization ⃗ ϵ 2 ⊥⃗ ϵ 1 , as specified in Eq. (6). Although defining a superradiance potential for the combined two-polarization problem is not straightforward due to the noncommutativity ofD AE 1 andD AE 2 , the singlemode potential of Eq. (36) can be defined for either of the two polarizations separately.
The definition of the potentials follows the same steps as in Sec. VI. First, we find a basis in which the dipole operators have the multi-two-level form of Eq. (24), i.e.,D þ 1 ¼ P α c ασeðαÞgðαÞ andD þ 2 ¼ P βcβσẽðβÞgðβÞ . The fgðαÞ; eðαÞg and fgðβÞ;ẽðβÞg bases are generally different, and they define two distinct sets of Bloch vectors ðS x α ; S y α ; S z α Þ and ðS x β ;S y β ;S z β Þ through Eq. (25), respectively. These Bloch vectors will evolve according to Eq. (29), and their dynamics can be parametrized by a single angular variable as in Eq. (32), S ⊥ α ðtÞ þ iS z α ðtÞ ¼ r α e i½c α θðtÞþφ α and S ⊥ β ðtÞ þ iS z β ðtÞ ¼r β e i½c βθ ðtÞþφ β . Notice that c α , r α , and φ α will generally be different fromc β ,r β , andφ β because they are associated with different two-level partitions of the multilevel structure.
When only one of the polarizations is present, the time evolution is described by [see Eq. (35)] ðd=dtÞθ ¼ −NΓðd=dθÞVðθÞ and ðd=dtÞθ ¼ −NΓðd=dθÞUðθÞ, respectively. Here, VðθÞ is defined as in Eq. (36), and analogously we have We call UðθÞ the orthogonal superradiance potential. For VðθÞ, we can set θ 0 ¼ jΩjτ as in Eq. (34) because of its relation to the Rabi preparation. However, for UðθÞ we set θð0Þ ¼ 0 for simplicity. The VðθÞ and UðθÞ potentials cannot be used to determine the time evolution of the two-polarization model (except when one polarization dominates). Nevertheless, their derivatives do reveal whether a given MF state is locally stable or unstable to decay withD − 1 andD − 2 , respectively.

B. Numerical results
We study the case of eight-level atoms with F g ¼ F e ¼ 3=2. Throughout this section, we work in the atomic basis defined with respect to the vertical quantization axisẑ V , and we decompose the decay operators into the verticalΠ AE and horizontalΣ AE components [see Fig. 2 and Eqs. (8) and (9)]. We consider two different initial ground states for the atoms: The atoms are excited with a Rabi drive of pulse area θ 0 and horizontal polarization ⃗ ϵ H , i.e., with theΣ AE operator. In thê z V basis [Eq. (9)], this operator connects the above ground states to je −1=2 i V and je 3=2 i V , creating a sort of zigzag motion in the level structure (see Fig. 10). Both of these excited states can decay via the horizontal polarizationΣ − , as well as through the vertical polarizationΠ − . In Figs. 10(a) and 10(e), we show the singlemode superradiance potential VðθÞ (red) associated witĥ Σ AE , which is the polarization of the Rabi drive. The potential for the two initial conditions jg 1=2 i V and jg −3=2 i V is markedly different. The potential for the initial state jg 1=2 i V [ Fig. 10

(a)] is given by VðθÞ
Þ and has two minima, i.e., two MF dark states. However, the potential for jg −3=2 i V [ Fig. 10(e)] is given by VðθÞ ¼ 1 2 ½1 − cosðθ= ffiffiffiffiffi 15 p Þ 3 and has instead saddle points, which, at the MF level, are semistable, i.e., stable to fluctuations in one direction (θ þ δ) but unstable in the opposite direction (θ − δ), 0 < δ ≪ 1. Note that in both cases the superradiance potential is periodic. This is due to the commensurability of the Clebsch-Gordan coefficients for linearly polarized transitions (C 0 m ∝ m), which leads toΣ þ ∝ simulations. Remarkably, for the jg 1=2 i V initial state [ Fig. 10(b)], all methods agree with each other, and the long-time value of the excited-state population matches the MF dark state (thin black dotted line) predicted by the single-mode potential VðθÞ. This suggests that the dynamics is indeed fully dominated byΣ − decay. This is confirmed in Fig. 10(d), which shows that light emission is happening almost exclusively in the ⃗ ϵ H polarization (light emission in the ⃗ ϵ V polarization is not exactly zero; see Sec. IX D). Moreover, Fig. 10(c) shows that the levels which were initially unpopulated, i.e., jg −1=2 i V , jg 3=2 i V , je −3=2 i V , and je 1=2 i V , remain approximately unpopulated to all times. While this is trivially true at the MF level, it seems counterintuitive when quantum fluctuations are included, since one would expectΠ − decay to populate such states.
When the atoms start in jg −3=2 i V , we obtain strikingly different excited-state population dynamics [ Fig. 10(f)]. For θ 0 ¼ 1.5π, all methods agree with each other, and the system decays back to the ground state jg −3=2 i V . However, for all other initial values of θ 0 shown, the MF dynamics does not match TWA and cumulant. In these latter cases, the MF evolution converges at long times toward the singlemode MF dark state of VðθÞ at n e ¼ 0.5. The beyond-MF dynamics, on the other hand, shows a two-step process. At early times, both TWA and cumulant appear to converge to the same single-mode dark state as MF; however, at longer times, quantum fluctuations make the system deviate from MF until it eventually settles on a different (quantum) dark state. When this happens, we observe decay in both theΣ − andΠ − decay channels, as demonstrated in the light emission profile of Fig. 10(h). This is further accompanied by a transfer of population into the initially unpopulated levels [ Fig. 10(g)]. Interestingly, the light emission shows various pulses of light being emitted at different times in different polarizations, signaling that decay in one polarization can trigger decay in the other polarization. All these results show that in this situation, both polarizations are playing an essential role in the dynamics.

C. Interpretation of results: Polarization competition
In order to understand the results of Fig. 10, we analyze the properties of the orthogonal potential UðθÞ associated withΠ − decay in this case. The key to understanding the findings is the realization that in MF approximation no decay with an operatorD − takes place if it has no initial coherences, i.e., if hD AE i ¼ 0. In our example, the initial pulse of ⃗ ϵ H light generates coherences in the drive's operatorΣ AE unless we end up at an extremum of VðθÞ. However, given our choice of initial ground state, we always have hΠ AE i ¼ 0 at t ¼ 0. This means that, regardless of the value of θ 0 , the initial state we prepare in our protocol is necessarily a MF dark state with respect toΠ − decay, i.e., ðd=dθÞUðθÞj˜θ ¼0 ¼ 0. The stability of these states to quantum fluctuations is then determined by the curvature of UðθÞ.
To evaluate UðθÞ, note first that the operatorΠ AE has the multi-two-level form of Eq. (24) when written in the V atomic basis as given in Eq.  37)], we can then compute the curvature of the orthogonal potential as a function of the initial pulse area θ 0 , which is the one that determinesr β . For jg 1=2 i V , we obtain ðd 2 =dθ 2 ÞUj˜θ ¼0 ¼ 1 120 Figures 10(a) and 10(e) show in blue the value of ðd 2 =dθ 2 ÞUj˜θ ¼0 obtained for each initial state θ ¼ θ 0 on the VðθÞ potential curve. Recall that if the curvature of UðθÞ is positive (negative), the state is stable (unstable) tô Π − decay. For the initial state jg 1=2 i V , we find that ðd 2 =dθ 2 ÞUj˜θ ¼0 > 0 in the vicinity of the two minima of VðθÞ. This proves that the MF dark states of VðθÞ are stable with respect to bothΣ − andΠ − decay.
At initial time, the values θ 0 ¼ ð1; 1.5; 3.5; 4Þπ chosen in Fig. 10(b) lie in regimes where ðd 2 =dθ 2 ÞUj˜θ ¼0 < 0. However, we do not observe a significant decay withΠ − at early times because of the superradiant delay t D , which we explain in Sec. VII. Initially, theΣ − channel is seeded (hΣ AE i ≠ 0), so the atoms can start decaying without delay NΓt D ∼ 1; cf. Eq. (38). However, since hΠ AE i ¼ 0 at t ¼ 0, coherences in this decay channel need to be built up, which leads to an NΓt D ∼ log N delay of theΠ − superradiant decay; cf. Eq. (39). Thus, in the limit of N → ∞, by the timeΠ AE coherences form, the atoms have long decayed viâ Σ − into one of the minima of VðθÞ, which are stable with respect to both polarizations. This is analogous to the competition between Clebsch-Gordan coefficients and initial coherences shown in Fig. 7.
For the initial state jg −3=2 i V , we find two distinct types of behavior. For values of θ where VðθÞ < 0.5, the state is stable toΠ − decay [ðd 2 =dθ 2 ÞUj˜θ ¼0 > 0]. This explains why the θ 0 ¼ 1.5π case is well described by MF and a single polarization. For values of θ where VðθÞ > 0.5, the state is unstable toΠ − decay [ðd 2 =dθ 2 ÞUj˜θ ¼0 < 0]. In these cases, theΣ AE decay initially dominates due to the initial coherences. However, because of the shape of the potential, the system always stays within the VðθÞ > 0.5 region, which is unstable toΠ − decay. Therefore, even though thê Π − decay is delayed by NΓt D ∼ log N, it will inevitably happen [as shown by the delayed superradiant pulse of light in Fig. 10(h)].
This explains the two-step decay process observed in Fig. 10(f): First, the system decays withΣ − toward the MF dark state of VðθÞ, and later it starts decaying withΠ − . This leads in Figs. 10(e)-10(h) to a complex dynamical interplay between the two polarizations, which eventually brings the system to a dark state that is not predicted by the single-mode VðθÞ potential. Moreover, the results of Figs. 10(a)-10(d) confirm the conclusion of Sec. VIII: Stable two-polarization MF dark states become asymptotically quantum dark states as N → ∞.

D. Light emission with orthogonal polarization
As we argue above, the collective decay is in some cases dominated by one of the polarizations. However, it is important to note that if we start in a product state, we always have light emission in both polarizations (unless the atoms are forbidden to emit in some polarization at the single-particle level, as for the six-level example of Sec. VIII). The reason for this is that we have both hΣ þΣ− i > 0 and hΠ þΠ− i > 0 for unentangled states; see Sec. IV D. Therefore, even if we have a MF dark state with respect to some polarization, light will be emitted at least with an intensity scaling with N, analogous to Sec. VIII B.
We demonstrate this in Fig. 11, which shows the light emission in both polarizations computed with cumulant for the jg 1=2 i V initial ground state. We show results for different atom numbers N with different markers. Figure 11(a) shows the results for an initial θ 0 ¼ 4π pulse, which corresponds to starting on the slope of VðθÞ; see Fig. 10(a). As demonstrated by the collapse of the curves when rescaling as hΣ þΣ− i=N 2 (red), the atoms emit superradiantly with an intensity scaling with N 2 in the excitation drive's polarization. As shown by the hΠ þΠ− i=N curves (blue), the atoms also emit into the orthogonal polarization, except with an intensity scaling with N. In comparison, Fig. 11(b) shows the light emission for a system starting at the MF dark state at θ 0 ¼ 2.542π. In this case, the light emission scales with N in both polarizations as the MF dark state decays into a quantum dark state; see Sec. VIII B.
While here we considered initial states with hΠ AE i ¼ 0, it is worth noting that one can easily prepare states with nonvanishing coherences in both polarizations. In such cases, we find that the dynamics is typically well described by MF for large N, usually shows superradiance scaling with N 2 in both polarizations, and often ends up in a dark steady state. Such two-polarization dark states can in principle be numerically found as described in Appendix I.

E. Delay time revisited
In Fig. 10, we saw an example of a superradiance potential VðθÞ with a saddle point. More generally, VðθÞ can contain extrema where arbitrary higher-order derivatives vanish. At these extrema, the dynamics induced by quantum fluctuations is extremely slow. To see this, let us assume that the potential has an unstable extremum at θ ¼ θ e and that the Taylor expansion around it gives to lowest nonvanishing order where n ∈ N. For θ ≪ 1, the dynamics is then described by _ θ ∝ NΓðθ − θ e Þ 1þn [cf. Eq. (35)].
In the two-level case and in most cases discussed above [see Figs. 8 and 10], the unstable extrema of the potential have n ¼ 0 because they correspond to maxima of simple sine and cosine functions. Solving the above equation for θð0Þ − θ e ≡ δθ 0 , jδθ 0 j ≪ 1 then gives θðtÞ ¼ θ e þ δθ 0 e NΓt . For a system starting exactly at the extremum, quantum fluctuations are of order δθ 0 ∼ ð1= ffiffiffiffi N p Þ. Solving for θðt D Þ ¼ const, we arrive at the well-known t D ∼ ð1=NΓÞ log N delay time of Eq. (39). In contrast, for n > 0 we obtain the solution θðtÞ ¼ θ e þ ½δθ −n 0 − nNΓt −1=n . This implies that the delay time scales polynomially as Thus, higher-order extrema lead to extremely slow dynamics at early times. Specifically, the saddle point of Fig. 10 is of order 2 þ n ¼ 3, which gives t D ∼ ð1=NΓÞ ffiffiffiffi N p .

X. EXPERIMENTAL IMPLEMENTATION
The multilevel dark states predicted here should be readily observable in a large variety of cavity QED settings. Alkaline-earth(-like) atoms trapped in optical cavities are particularly well suited for two reasons. First, alkaline-earth atoms have a relatively simple electronic structure which contains a unique ground-state manifold. This simplifies the physics, as it reduces the number of states to which an excited state can decay.
Second, alkaline-earth atoms feature ultranarrow lowlying electronic transitions [137] between 1 S 0 and either 3 P 0 or 3 P 1 . For example, the 1 S 0 ðF ¼ 1=2Þ to 3 P 1 ðF ¼ 3=2Þ transition in 171 Yb is equivalent to the six-level model considered above, whereas the 1 S 0 ðF ¼ 9=2Þ to 3 P 0 or 3 P 1 ðF ¼ 9=2Þ transition in 87 Sr is a 20-level generalization of the eight-level system presented above. When the cavity is coupled to one of these narrow transitions, both the collective dynamics induced by superradiance and the timescales at which dark states are stable (before they decay due to single-particle emission; see Sec. XI) cover a range of timescales that can be experimentally accessed in standard cavity QED settings.
Alternative implementations can be done in arrays of alkali atoms featuring dipole-allowed transition via Raman dressing. In this case, the cavity is effectively coupled to an optically dressed atomic ground-state manifold that is engineered to imitate a long-lived optically excited atomic state similar to alkaline-earth atoms. This protocol was used in Ref. [74] to engineer an effective two-level system but FIG. 11. Light in orthogonal polarization. Cumulant results are shown for the eight-level system with F g ¼ F e ¼ 3=2 starting in jg 1=2 i V and being initially excited by aΣ AE drive. We plot the light intensity in both polarizations, hΣ þΣ− i (red) and hΠ þΠ− i (blue), for different atom numbers N (markers shown in legend). The different polarizations are scaled differently with N as indicated in the plots. The collapse of the curves demonstrates the corresponding N scaling. Panel (a) shows results for θ 0 ¼ 4π, which is a nondark state. Panel (b) shows θ 0 ¼ 2.542π, which corresponds to a minimum of VðθÞ, i.e., a MF dark state; see potential in Fig. 10(a).
can be generalized to emulate the multilevel configurations discussed here.
In these systems, the protocol proposed to observe superradiance and dark states relies simply on (i) preparing the atoms in a well-defined ground state, and (ii) exciting them with a pulse of the desired polarization, which can be easily implemented. In particular, optical pumping is a simple way to prepare all atoms in the same ground state, typically in the stretched states jg AEF g i. Experiments can also straightforwardly rotate this ground state by using microwave lasers or magnetic fields.
As an example, we consider the 20-level system F g ¼ F e ¼ 9=2, which is relevant for 87 Sr. We consider the same scheme as we study for the eight-level system in Sec. IX. We assume that the atoms start in the jg −9=2 i V ground state, which can be easily prepared via optical pumping, and are then excited with horizontal polarization Σ AE . This leads to a particularly simple looking superradiance potential VðθÞ Fig. 12(a), which looks quite similar to the eight-level potential of Fig. 10(e). In particular, we find that there is also a saddle point but with an extremely flat potential around it. The same is true for the curvature of the orthogonal potential, which is given by ÞÞ. Similar to the eight-level case, the saddle point is overall unstable, either toΣ − orΠ − decay. However, due to the flatness of the potentials, the states close to the saddle point are metastable for a very long time before they decay to their true final state. Specifically, the superradiance delay time scales as t D ∼ ð1=NΓÞN 7=2 according to Eq. (46). Because of this, the saddle-point states will appear dark to cavity decay for experimentally relevant timescales.
The slow time evolution around the saddle point is illustrated in Fig. 12(b). The closer to the saddle point, the longer it takes for the states to start decaying. For initial states far away from the saddle point, the decay is analogous to the eight-level case of Sec. IX. For states with Vðθ 0 Þ < 1=2, the atoms decay back to the initial ground state, jg −9=2 i V . For states with Vðθ 0 Þ > 1=2, the atoms decay with both polarizations into a dark state that is not predicted by the single-polarization potential. The timescales of these processes scale as NΓ and should be experimentally observable. [Note that when the atoms start exactly at the global maximum of VðθÞ (θ 0 =π ¼ 3 ffiffiffiffiffi 11 p ), the system decays back to the ground manifold. This is because the initial state is je 9=2 i V and the Clebsch-Gordan coefficient for e 9=2 → g 9=2 is larger than for e 9=2 → g 7=2 , so the atoms will predominantly decay viaΠ − into g 9=2 ; cf. Sec. VII. The same reasoning applies to the eight-level example of Fig. 10 (lower panel) when the atoms start at the global maximum of VðθÞ, which corresponds to the je 3=2 i V state.]

XI. ROBUSTNESS OF DARK STATES IN EXPERIMENTS
We anticipate three possible sources of experimental imperfections, which were not included in our model so far: (i) stray magnetic fields, (ii) inhomogeneous couplings, (iii) coherent cavity interactions, and (iv) spontaneous decay into free space. We analyze in this section whether dark states are observable when these additional processes are included.

A. Magnetic fields
First, the magnetic field in experiments is never exactly zero. Nonvanishing magnetic fields can lead to mixing between the internal levels and, thus, to decay of the dark state. However, if the induced Zeeman splittings δ g;e of the g and e manifolds are small compared to the superradiant decay rate ΓN, dark states can live for a long time and even modified dark states can emerge. This means that the effect of δ g;e can be suppressed by using a large number of atoms N.
To show this, we investigate in Fig. 13(a) the 20-level example of Sec. X (see Fig. 12) including now Zeeman shifts through the HamiltonianĤ B ¼ δ g P m mσ V g m g m þ δ e P m mσ V e m e m with δ g ¼ 2 3 δ e ≡ δ z . As an example, we consider the case θ 0 ≈ 5π where the atoms start in the saddle point of VðθÞ [ Fig. 12(a)]. For δ g;e ¼ 0, the atoms simply stay in this metastable dark state for all times shown due to the long delay time, t D ∼ ð1=NΓÞN 7=2 . For a small Zeeman shift δ g;e ≪ ΓN, the system stays in this metastable state for some time but then slowly starts decaying at a reduced decay rate due to the Zeno effect [138]. Notice that the smaller δ g;e , the longer it takes for the system to start decaying. After the initial decay, the system eventually settles down into a new dark state. For large enough magnetic fields (δ g;e ∼ ΓN), we find that the system superradiantly decays to the ground manifold.

B. Inhomogeneous couplings
In typical experiments, the atoms are tightly trapped by a deep one-dimensional optical lattice that is supported by an optical cavity. Because of the wavelength mismatch between the lattice laser and the cavity field, the coupling of the light to the atoms varies between lattice sites. This setup can be modeled by a modified master equation with where the operatorsD AE i;γ are now collective operators as defined in Eq. (2) but acting only at lattice site i. The inhomogeneity of the interactions enters through the coefficients ξ j ¼ cosðwjÞ with w ¼ πλ L =λ c , where λ L is the lattice laser wavelength and λ c the cavity mode wavelength [32].
Despite the inhomogeneity, we can perform an analogous single-mode analysis as in Sec. VI. For this purpose, we define for each of the sites i an independent set of α-Bloch spheres described by the Bloch vectors ⃗ S i;α . All these Bloch vectors couple to each other through a modified total dipole moment ⃗ D → P i ξ i ⃗ D i , where ⃗ D i is the dipole moment at site i. Following the same derivation as in Sec. VI, one can show that collective decay is then determined by a superradiance potential that combines the dipoles of all these two-level systems as FIG. 13. Magnetic field and inhomogeneous couplings. Results for the 20-level system with F g ¼ F e ¼ 9=2 starting in the jg −9=2 i V state and being initially excited by aΣ AE drive of pulse area θ 0 [cf. Fig. 12(a)]. (a) Time evolution of the excited-state population n e ¼ 1 N P m hσ q e m e m i for an initial pulse θ 0 ¼ 5π lying at the saddle point of VðθÞ. Shown are MF (solid) and TWA (dashed) simulations averaged over 10 4 trajectories for N ¼ 10 4 and for different Zeeman splittings δ z =ðNΓÞ with δ g ¼ 2 3 δ e ≡ δ z (motivated by Ref. [139]). Note that dashed lines are not visible because they lie on top of the solid MF lines of the same color. For small enough δ z =ðNΓÞ, the system decays into a dark state. (b) Superradiance potential VðθÞ forΣ − decay (upper panel) and curvature of the orthogonal potential ðd 2 =dθ 2 ÞUðθÞj˜θ ¼0 forΠ − decay (lower panel) for a system with inhomogeneous couplings. We show results for different number of lattice sites (color lines) computed for λ L ¼ 813 nm and λ c ¼ 689 nm [32]. The one-site case (i ¼ 0) corresponds to Fig. 12(a). For more sites, the potential is modified but still contains stable dark states.
Here, we assume that atoms in all lattice sites start in a state with the same initial r α , φ α . The sum over lattice sites modifies the shape of the superradiance potential VðθÞ computed in previous sections. However, the minima of this modified potential VðθÞ will still correspond to (single-mode) MF dark states of the dynamics. These dark states are now states where the sum of the dipole moments for each internal transition and lattice site cancels out. Figure 13(b) shows how the superradiance potential for the F g ¼ F e ¼ 9=2 example of the previous section is modified as we include more lattice sites into the system. Although VðθÞ becomes more complicated due to the higher number of frequencies, it still contains a large number of stable MF dark states.

C. Coherent interactions χ
As we argue in Sec. IV, the coherent part of the cavitymediated interactions χ can be set to zero by working on resonance, Δ c ¼ 0 [Eq. (7)]. When χ ≠ 0, dark states still remain dark, as shown in Sec. IV; however, the superradiant dynamics is significantly altered. As an example, we show in Fig. 14 the dynamics for the 20-level example of Fig. 12 for different values of θ 0 , with χ ¼ Γ (circles) compared to χ ¼ 0 (solid line). When the system starts at the (extremely flat) saddle point of the potential θ 0 ≈ 5π, it remains in this quasidark state with and without χ for all times shown. In contrast, for values of θ 0 on the slope of the potential VðθÞ, the dynamics is considerably modified by χ. Nevertheless, the system still ends up in a dark state, although a different one. Note, in particular, that for θ 0 ¼ 2π the system goes back to the ground state unless χ ≠ 0.
This example shows that to describe the dark states and dynamics in the presence of χ, the simple picture offered by the superradiance potential VðθÞ would need to be extended. The rich and intricate dynamics induced by χ will be explored in subsequent work.

D. Spontaneous decay and dipolar interactions
A final important source of decay for the dark state is single-particle spontaneous decay γ s into free space. The dark states presented here are only dark with respect to decay into the cavity modes. Hence, spontaneous decay will inevitably lead to the decay of such dark states. The relevant parameter is the ratio γ s =ΓN with respect to the collective decay rate. Thus, increasing ΓN is one way to make γ s negligible for the relevant timescales.
Similarly, vacuum-mediated dipolar interactions between the atoms can break the collective nature of the system and will generally couple dark states to states that radiate. Dipolar interactions can be large if the atoms are placed at a close distance r since they scale as γ s =ðk 0 rÞ 3 at short distances, where k 0 ¼ ω 0 =c. However, the typical interparticle distances in an optical cavity experiment are of order of the transition wavelength r ∼ λ 0 =2 (λ 0 ¼ 2π=k 0 ) or larger. In this regime, k 0 r ≳ 1, so dipolar interactions will be of order γ s and will also be suppressed compared to ΓN.

XII. CONCLUSION AND OUTLOOK
The main result of this work is that, even within the permutationally symmetric (collective Dicke) manifold, superradiant dynamics in multilevel atoms can get stuck in dark states that retain a finite fraction of excitations and are entangled (see Sec. I A for a summary). These dark states open up a number of exciting research directions.
One of the potential applications of dark states is in quantum sensing and metrology. Dark entangled states within the collective manifold will not suffer from superradiant decay and will therefore retain for longer time the phase information necessary for Ramsey-type protocols. An interesting question in this regard is how correlations will evolve when the system starts close to or at a meanfield dark state, especially in the presence of nonzero coherent interactions χ. This could lead to the creation of multilevel entangled squeezed states [107,140].
The complexity of multilevel systems, even when constrained to a collective manifold, also anticipates a rich landscape of physical phenomena. For example, dark states can lead to interesting phases of driven-dissipative systems, since they can give rise to multiple steady states. Apart from this, long-lived subradiant states can also be useful for storage of quantum information and for the implementation of quantum repeaters [141]. Moreover, the complexity of the collective eigenstate structure could be useful for engineering nontrivial effective FIG. 14. Dynamics in the presence of elastic-cavity-mediated interactions, χ ≠ 0. Results for the 20-level system with F g ¼ F e ¼ 9=2 starting in the jg −9=2 i V state and being initially excited by aΣ AE drive [cf. Fig. 12]. Time evolution of the excited-state population n e ¼ ð1=NÞ P m hσ q e m e m i for different initial pulse areas θ 0 . We compare results for χ ¼ 0 (solid) with χ ¼ Γ (circles) computed with TWA for N ¼ 10 4 and 10 4 trajectories. We checked that the results agree with MF. When starting at the very flat saddle point of the potential (θ 0 ≈ 5π), χ has no effect on the population dynamics. In the other cases, χ makes the system decay into a different dark state.
Hamiltonians in the ground-state manifold analogous to Raman dressing setups [142].
Finally, one particularly appealing ingredient is that the cavity-mediated interactions can be arranged to match characteristic collisional timescales. Investigating the interplay of both contact and infinite-range cavity-mediated interactions [143], especially within the long-lived darkstate manifolds investigated here, could open untapped opportunities for the simulation of rich models of orbital quantum magnetism, novel phases of matter, chaotic dynamics, and for the engineering of minimal models of holographic gravity [39,40].

APPENDIX A: BASIS ROTATIONS
In this section, we provide details on how to transform between the different atomic bases introduced in Sec. II and Fig. 2. In general, we consider angular momentum eigenstates jj; mi q ≡ jmi q defined with respect to a quantization axisẑ q . We can transform between different quantization axes by applying a unitary rotation jmi q 2 ¼Rjmi q 1 (R −1 ¼R † ) which rotatesẑ q 1 →ẑ q 2 . Note that for a generic state jψi ¼ P m ψ q m jmi q , the coefficients transform instead as ψ q 2 m ¼ R † mn ψ q 1 n with R mn ¼ hmj qR jni q (repeated indices are summed over). Similarly, forÔ ¼ P m;n O q mn jmihnj q we obtain O q 2 mn ¼ R † mk O q 1 kl R ln . Adopting an active convention [144], a generic rotation operator is given byRðφ; θ; χÞ ¼ e −iφĴ q z e −iθĴ q y e −iχĴ q z , where J q x;y;z are spin-j generators fulfilling the usual SU(2) commutation relations and are also defined with respect to some axisẑ q . Given the quantization axis reference frames defined in Fig. 2, i.e., fx H ;ŷ H ;ẑ H g ¼ fx V ; −ẑ V ;ŷ V g and fx k ;ŷ k ;ẑ k g ¼ fẑ V ;ŷ V ; −x V g, the corresponding atomic bases can then be transformed as Notice that when performing these rotations, the states acquire global phase factors, which in some cases can be dropped for simplicity.
The multilevel atoms considered in this work have a level structure composed of a ground-and an excited-state manifold with angular momenta F g and F e , respectively. The transformations of Eqs. (A1) and (A2) can be easily generalized to this case by considering the rotation generators in the spinðF g Þ ⊕ spinðF e Þ representation, i.e., where α ¼ x, y, z. Written as a matrix, J q α is a block diagonal matrix where the first block isĴ qðF g Þ α and the second oneĴ

APPENDIX B: NUMERICAL METHODS
In this section, we provide some further details on the different numerical approximations used in the main text, which are introduced in Sec. III. For simplicity, we drop in this section the label q that denotes the quantization axis employed.

Mean-field approximation
In the limit of large N, quantum fluctuations and correlations of collective systems are typically suppressed by 1=N. A common approach is then to use a MF approximation which neglects two-body correlations. The equations of motion for the single-body expectation values hσ ab i depend on two-body expectation values hσ abσcd i, whereσ ab andσ cd and a, b, c, d stand for any internal level of the atoms. In order to obtain a closed set of MF equations for the single-body observables, we approximate This results in approximately l 2 MF equations of motion for the one-point expectation values hσ ab i, which can be efficiently solved numerically. The equations obtained from using Eq. (B1) are the ones employed in this paper. However, it is important to note that this decoupling scheme is not unique. A common point of ambiguity in the MF decoupling of products of collective operators hσ abσcd i ¼ P i;j hσ i abσ j cd i is how to treat the "self-interaction" i ¼ j terms, where i, j label a single atom. In the approximation of Eq. (B1), we effectively assume hσ i abσ i cd i ≈ hσ i ab ihσ i cd i. Another alternative, however, would be to simplify the product at the quantum level as hσ i abσ i cd i ¼ δ bc hσ i ad i, where δ bc is a Kroenecker symbol.
This issue has been discussed in related works [132,145,146] for two-level systems. The conclusions seem to depend on the problem at hand, but generally speaking, the method of Eq. (B1) appears to be more accurate when the system remains in the permutationally symmetric manifold of pure states, i.e., when processes such as single-particle spontaneous emission are absent. At large N, we checked that the differences in the time evolution between the two schemes are negligible (order 1=N). However, in our superradiant problem, we find that using the second decoupling scheme hσ i abσ i cd i ¼ δ bc hσ i ad i leads to an unphysical decay of the dark states with a rate of order 1=N, which is absent in ED simulations (we note that this behavior is consistent with the idea that these selfinteraction terms act as effective single-particle decay terms in the hσ i abσ i cd i ¼ δ bc hσ i ad i decoupling scheme [132]). Because of this, we employ Eq. (B1) in both MF and TWA simulations.

Truncated Wigner approximation
One simple way to take quantum fluctuations into account is by employing TWA. TWA consists in solving the MF equations but with the initial conditions sampled from a noise distribution that models the initial quantum fluctuations of the system [147]. Expectation values are then obtained by averaging over MF trajectories. TWA usually works when the quantum-noise distribution does not spread too much over the phase space, but it has not been well tested for dissipative system so far. Nevertheless, it can provide us valuable insights into the role of quantum fluctuations.
For the TWA simulations, we employ the MF equations obtained from Eq. (B1). Typically, the TWA equations of motion should be obtained using instead a symmetric decoupling scheme [147] 1 2 hfσ ab ;σ cd gi ≈ hσ ab ihσ cd i. However, the difference between the symmetric scheme and Eq. (B1) lies only in self-interaction terms of the sort described in the previous subsection, which we neglect for the same reasons as for MF.
For the sampling over initial conditions, we employ a multivariate Gaussian distribution as detailed in Ref. [148]. This continuous approximation reproduces the mean and variance of the quantum-noise distribution and is justified in the limit of large N, where it should be equivalent to discrete sampling schemes [149]. Specifically, we define a complete set of single-body Hermitian operators composed of the Gell-Mann-type operatorsσ x ab ¼ 1 2 ðσ ab þσ ba Þ and σ y ab ¼ 1 2i ðσ ab −σ ba Þ for a > b, andσ aa for a ¼ b. We denote these operators asÔ k , where k runs from 1 to l 2 . For each operator, we compute the quantum expectation values μ k ≡ hÔ k i and C kq ≡ 1 2 hfÔ k ;Ô q gi for all k, q. For each TWA trajectory, we initialize the classical variables O k associated with theÔ k operators by drawing a set of random numbers from a multivariate Gaussian distribution with mean μ k and covariance matrix Σ kq ≡ C kq − μ k μ q . This can then be transformed into any other singlebody basis.
It is important to note that when computing C kq from two-body quantum expectation values, the self-interaction terms must be treated exactly; i.e., we use hσ i abσ i cd i ¼ δ bc hσ i ad i. Even though such terms are subleading in 1=N, they are essential in cases where the mean dipole is zero such as when starting at an unstable MF dark state.

Cumulant expansion
Another way to include the effect of correlations and quantum fluctuations is by employing a cumulant expansion which takes higher-order correlators into account. Specifically, we include both one-body hσ ab i and two-body correlators hσ abσcd i and derive equations of motion for these quantities by neglecting the connected part of threebody correlators as The number of cumulant equations of motion obtained in this way scales as l 4 . This approximation should work as long as connected three-point correlations are small, and it is used to compare to TWA predictions.
As in the mean-field case, the decoupling scheme of Eq. (B2) is not unique. If we treat self-interaction terms such as hσ i abσ i cdσ j ef i ¼ δ bc hσ i adσ j ef i exactly, we obtain additional terms to Eq. (B2) which are suppressed by 1=N and 1=N 2 compared to the rest. However, we checked that such terms lead again to a decay of the dark state with a rate of order 1=N 2 , which sometimes settles in a different dark state. Because of this, we employ instead Eq. (B2).

APPENDIX C: EIGENSTATES
In this section, we complement the eigenstate analysis of Sec. IV with a discussion of the conservation laws of the system and the eigenstates of level structures not covered in the main text. In particular, we show that the size of the eigenvalue problem for finding eigenstates with k excitations can be reduced to diagonalizing d × d matrices with at most d ∼ k 2F e N maxð0;2F g −2Þ . These simplifications can help push numerical simulations of future studies of collective multilevel systems beyond the naive d ∼ N l−1 scaling of the collective Hilbert space; see Eq. (14).

Conserved quantities
When we write the effective Hamiltonian Eq. (15) in thê z k basis, we can straightforwardly identify several interesting conservation laws which would have been otherwise hidden in a different basis.
First, if we connect the atomic states in theẑ k basis in a zigzag motion [see, e.g., Fig. 10], we can define two sets of levels (A, B) which conserve particle number independently from each other. To formalize this, we define the magnetic number m as even (uneven) if the integer m þ F g is even (uneven). The conserved number operators for the A, B sets can then be defined as whereN A þN B ¼ N and we defineN q g m ¼σ q g m g m , N q e m ¼σ q e m e m . [We note that another way to identify the A, B sets of states is by defining a parity operatorP ¼ e iπF g e iπĴ k z (see Appendix A). All even (uneven) jg m i k and je m i k states are then þ1 (−1) eigenstates ofP. It can be straightforwardly shown that fP;L AE g ¼ fP;R AE g ¼ 0; i.e., L AE andR AE connect states of opposite parity. This justifies the definitions of Eqs. (C1)  Another conserved quantity is the total magnetic number in theẑ k atomic basis defined bŷ This can be easily seen as bothL þL− andR þR− separately conserveM k . Apart from this,Ĥ eff has a discrete m → −m symmetry, where m is the magnetic number m in theẑ k basis. The latter symmetry is because of the jC þ1 m j ¼ jC −1 −m j symmetry of the Clebsch-Gordan coefficients.
In Sec. III, we show that the size of the permutationally symmetric Hilbert space scales roughly as N l−1 with l ¼ 2ðF g þ F e þ 1Þ. Using the conservation of the total excitationsN e (see Sec. IV), the particle numberN A (or N B ), and the magnetizationM k , we can thus reduce the eigenvalue problem to diagonalizing d × d matrices where d scales at most as N l−4 .

Four-level systems
In Sec. IV B, we argue that the eigenstates of four-level systems are simply given by the PS states of Eq. (12) in thê z k basis. For the F g ¼ 0, F e ¼ 1 level structure, the PS states can be parametrized as where N A þ N B ¼ N, N A ≥ k L þ k R . The dipole operators in this case are given byL − ¼σ k g 0 e −1 andR − ¼σ k g 0 e 1 . Hence, the decay rates of the above PS states read Note that the je 0 i k state is decoupled from the dynamics at the single-particle level, and hence, all states with k L ¼ k R ¼ 0 are trivially dark. For the F g ¼ 1, F e ¼ 0 level structure, the PS states can instead be parametrized as where N B þ k þ n R þ n L ¼ N. In this case, we haveL − ¼ ð1= Here, jg 0 i k trivially decouples from the dynamics, but there are no single-particle dark states. The decay rates of both Eqs. (C5) and (C7) are in essence a sum over two two-level decay rates; cf. Eq. (16). Note that in the two level structures above, theL AE and R AE do not commute (in contrast to F g ¼ F e ¼ 1=2), but thê L þL− andR þR− do.

Six-level systems
For systems with six levels (or more), theL AE andR AE operators involve at least two different transitions each (see Sec. IV C). Because of this, the productsL þL− andR þR− appearing in the effective Hamiltonian [Eq. (15)] have cross terms that are responsible for PS states [Eq. (12)] not being eigenstates. Instead, acting withĤ eff on a PS state results in a superposition of other PS states weighted with combinations of Clebsch-Gordan coefficients. Because of this, some eigenstates can become dark to cavity decay, as we discuss in the main text.
In Fig. 15, we show again histograms of the decay rates of eigenstates for fixed excitation number k but this time for different six-level systems with (red) Alongside, we also provide the percentage of eigenstates that are dark for each k and level structure (in the corresponding color). As one would expect, the plots show that the fraction of dark states is larger the more excited levels there are, i.e., the fewer decay channels there are. For analogous reasons, the fraction of dark states decreases with k. The Λ-type ðF g ; F e Þ ¼ ð3=2; 1=2Þ level structure is a special case where we find there are no dark states.
Despite the added complexity, the problem still has some structure left due to the symmetries discussed above. When written in the basis of PS states in theẑ k basis,Ĥ eff splits into blocks of size d × d with d scaling at most as N l−4 but often considerably better than that. For simplicity, we write In the following, we discuss the eigenstates of different sixlevel systems and elaborate on the scaling of the block sizes ofĤ ¼L þL− þR þR− with some examples. We note that, while we focus on states with a few number of excitations k, analogous arguments can be applied to states with N − k excitations. Note also that the bounds we discuss are not tight and are only intended as an order-of-magnitude estimate.
ffiffi ffi 3 p Þσ k g −1=2 e 1=2 . This means that and similar forL þL− . The terms in the second line of the above equation connect different PS states with each other.
In general, applying the operatorR þR− repeatedly to a state of the form j a b f c g d i k connects it to states of the form and analogously for L þL− . This leads to the rough N 2 ¼ N l−4 scaling of the blocks ofĤ. However, as we see above, for a small fixed number of excitations k ¼ a þ b þ c þ d, the size of the problem is independent of N. For fixed k, the size of the block is upper bounded by ða þ b þ 1Þðc þ d þ 1Þ ≲ ðk=2 þ 1Þ 2 ∼ k 2 , which is independent of N. To see this, note that repeated application of the cross terms ofR þR− is equivalent to a transfer of population e 3=2 → e 1=2 and g −1=2 → g 1=2 (and vice versa). This process can happen only while the populations are non-negative, which gives at most a þ b þ 1 distinct states, and analogously forL þL− .
Þðσ g 0 e 1 þσ g −1 e 0 Þ. As an example, we consider again a state with a single excitation, e.g., j 1 up to a normalization. In general, the eigenvalue problem for this level structure turns out to be of similar complexity as the previous one: For a fixed number of excitations k, the size of the blocks to be diagonalized is independent of N. For example, k ¼ 2 states such as j 2 a 0 b 0 c i k get mixed with five other states only. For fixed but general k, repeated application ofĤ can transfer population between any two excited states, as long as the ground-state populations allow it. Specifically, L þL− transfers e −1 ↔ e 0 andR þR− transfers e 1 ↔ e 0 . Therefore, the number of states that are connected byĤ is FIG. 15. Six-level collective eigenstates. Histograms of the decay rates (γ=Γ) of all collective eigenstates with (a) k ¼ 1, We provide in the plots the percentage of eigenstates that are dark for each k and level structure in the associated color. Note that the histograms are stacked on top of each other. upper bounded by how many ways there are to distribute k excitations among the three excited states. This is given by the binomial ð kþ2 2 Þ, which scales as k 2 .
c. F g = 3=2, F e = 1=2 For F g ¼ 3=2, F e ¼ 1=2, the single-excitation eigenstates cannot be computed independently of N as before. To see this, note first that the dipole operators are The crucial difference from the previous case is that now, even for a state with a single excitation, repeated application ofĤ can reshuffle the population of the ground states while keeping the excited-state population untouched. For example, applyingL þL−RþR− on j a We call this type of process where the ground population changes but the excited population does not (or vice versa) a loop process. Through repeated application of the loop process, the state j a 1 b 0 c d i k will connect to states j aþx dÞ ≤ x ≤ minðb; cÞ. The number of possible values of x is then minðb; cÞ þ minða; dÞ þ 1, which scales at most as N=2.
For a general number of excitations k, repeated application of onlyL þL− (orR þR− ) leads to a transfer of excitations e −1=2 ↔ e 1=2 , which generates ð kþ1 1 Þ ∼ k states with different population distributions of the excited states. Moreover, for each distinct excited-state distribution, repeated application of the loop process generates at most order N other states with different ground-state populations. Therefore, the size of the largest block inĤ scales as kN. Note that this rough scaling applies best for small k and does not take possible constant prefactors into account. In particular, the scaling can be shown to become independent of N when k ≈ N.

l ≥ 8-level systems
For systems with l ≥ 8 levels, the eigenvalue problem becomes increasingly complicated. Numerically, we find that dark eigenstates exist for any system with l ≥ 8 levels. In Fig. 16, we show histograms of the decay rates of eigenstates for different eight-level systems with (red) We show results for fixed excitation number k and N ¼ 6, and we provide the percentage of eigenstates that are dark for each k and level structure (in the corresponding color). As before, the fraction of dark states is larger for small k and for larger number of excited states. As opposed to the six-level case, however, here all level structures contain dark states.
To estimate the size of the eigenvalue problem, i.e., the size of the blocks ofĤ [Eq. (C8)] when written in theẑ k basis, we can proceed as in the previous section. Through repeated application ofL þL− andR þR− on a PS state, we can reshuffle the population of the excited states. In the worst-case scenario, all possible combinations of excitedstate populations with k total excitations are connected to each other. The total number of excited-state combinations scales at most as ð kþ2F e 2F e Þ ∼ k 2F e . Note that this scaling is more favorable for some l ≤ 6 cases, as we discussed in previous sections.
Each state with a given distribution of excited-state populations can be connected through repeated application ofĤ to other states with the same excited but different ground distribution. We call such processes loop processes, as we discussed in the previous section. Each loop process can be written as a combination of irreducible loop processes present inL þL−RþR− (andR þR−LþL− ). The number of independent irreducible loop processes, i.e., loops that cannot be written in terms of other loops, is maxð0; 2F g − 2Þ. To see this, note that, in order to leave the excited-state population untouched, each loop process in L þL−RþR− must involve transferring population from a pair of ground states ðg m ; g m 0 Þ to another pair ðg mþ2 ; g m 0 −2 Þ with m ≠ m 0 . Note that this is only possible if F g ≥ 3=2. If we fix m ¼ −F g , then each m 0 ≥ −F g þ 3 constitutes an FIG. 16. Eight-level collective eigenstates. Histograms of the decay rates (γ=Γ) of all collective eigenstates with (a) k ¼ 1, We provide in the plots the percentage of eigenstates that are dark for each k and level structure in the associated color. Note that the histograms are stacked on top of each other. independent loop process. All other loop processes with m; m 0 ≠ −F g can be obtained by combining m ¼ −F g loop processes.
Each independent loop process connects a given state with at most order N distinct states. Therefore, the size of the largest block ofĤ scales at most as k 2F e N maxð0;2F g −2Þ .
For F g ¼ 1, F e ¼ 2, we obtain a k 4 scaling which is again independent of N. For F g ¼ F e ¼ 3=2 we get k 3 N and for F g ¼ 2, F e ¼ 1 we get k 2 N 2 .

APPENDIX D: PROOF THAT QUANTUM DARK STATES MUST BE ENTANGLED
In this section, we present a proof showing that if the system is in a product state, and the internal level structure does not admit single-atom dark states, then the state cannot be a quantum dark state. This implies that many-body quantum dark states are necessarily entangled and, therefore, that the presence of nondecaying excitations is an entanglement witness.
To show this, we assume without loss of generality that the system can decay only via one polarizationD AE ¼ P iD AE i , where i labels the atoms. We assume that the system is in an arbitrary (possibly mixed) product state of the form where p k > 0 andρ ðkÞ i are single-atom density matrices. Without loss of generality, we can assume that ρ Importantly, we assume that at least one of the singleparticle jψ ðkÞ i i states fulfillsD − i jψ ðkÞ i i ≠ 0; i.e., at least one of the atoms in the mixture is excited and not in a single-particle dark state.
In the main text (Sec. VI) we show that a MF dark state must have hD AE i ¼ 0. A product state can fulfill this because i . This expression can be zero even if all atoms have nonzero dipole moment since the various terms in the sum P k can cancel each other out. Thus, we find again that a mixed product state can be a MF dark state. However, in the full quantum theory, a dark state also must have zero quantum fluctuations in the dipole moment operator. In particular, a quantum dark state has to fulfill hD þD− i ¼ 0. To see why, notice that the equation of motion for the total inversion S z ¼ 1 2 ð P mσemem − P nσgngn Þ is given by ðd=dtÞhŜ z i ∝ hD þD− i because ½D AE ;Ŝ z ¼∓D AE . Thus, if hD þD− i ≠ 0, then the total number of excitations is not stationary.
We show in the following that if the system is in a product state given by Eq. (D1) and it contains a nonzero amount of excitations (which are not single-particle dark, as specified above), then hD þD− i > 0. For this, we write The last inequality follows from assuming that at least one of the i, k terms in the sum fulfillsD In Sec. VI, we make a variable transformation from ðS x α ; S y α ; S z α Þ into ðS k α ; S ⊥ α ; S z α Þ which leads to the ansatz of Eq. (32) and the equation of motion of Eq. (33). The explicit transformation used rotates the α-spins such that one of the components is parallel to the torque vector ⃗ Ω as where ⃗ Ω ¼ ðΩ x ; Ω y ; 0Þ and ⃗ Ω ⊥ ¼ ð−Ω y ; Ω x ; 0Þ. This transformation can also be written as The solution to the Rabi equations is then given by Thus, the total population imbalance evolves as P α S z α ¼ P α r α sinðc α jΩjτ þ φ α Þ, which by comparison with Eq. (36) proves the relationship of Rabi oscillations to the superradiance potential, Eq. (37).
For arbitrary initial conditions that are not necessarily prepared by starting in a ground state and applying a Rabi drive, we can describe the decay dynamics by applying instead the transformation where ⃗ D ⊥ ¼ ð−D y ; D x ; 0Þ, ⃗ D ¼ ðD x ; D y ; 0Þ, and we generally set σ 0 ¼ 1. Note that for the excitation-decay scenario considered in the main text, the definitions of Eqs. (E1) and (E3) are identical if we set σ 0 ≡ sign½ P α c α S ⊥ α ð0Þ with S ⊥ α ð0Þ as defined in Eq. (E1).

APPENDIX F: DIPOLE OPERATOR FOR ELLIPTICAL POLARIZATION
In Sec. VI, we explain how the dipole operatorD AE can be brought into the multi-two-level form of Eq. (24) by choosing the appropriate quantization axis. In some cases, however, that strategy does not suffice. For the general case of elliptical polarization (perpendicular to the cavity axis), we can instead find the basis that diagonalizes the commutator ½D þ ;D − . The logic behind this is that if there exists a basis such thatD AE has the form of Eq. (24), then the commutator ½D þ ;D − has to be diagonal. [From a Lie algebra perspective, theD AE operators can be decomposed as a sum of root operatorsÊ ð †Þ α . The basis whereD AE can be written as in Eq. (24) is then equivalent to the basis where eachÊ ð †Þ α can be written as jgðαÞiheðαÞj (or jeðαÞihgðαÞj).] As an example, we consider the elliptical polarization obtained with ED simulations for N ¼ 30. The distribution is maximal for Δn g ¼ −1 and decays exponentially for increasing Δn g . The reason for this is that decaying into one ground state increases the probability of decaying the next time into the same ground state. For comparison, we also show in gray the binomial distribution one would obtain if the atoms decayed independently into the cavity with the decay probabilities of Eq. (G1).

Three-level system: Balanced case
An interesting variation of the previous example is to consider a three-level system with equal decay probabilities in both channels, This is the case if we consider the jg AE1=2 i V and je −1=2 i V states of the ðF g ; F e Þ ¼ ð1=2; 1=2Þ level structure. Figure 17(b) shows the resulting probability distribution for the ground-state population. While for independent atoms (gray) we get a binomial distribution centered on zero, for superradiant decay we obtain a flat distribution due to the collectively increased probability of decaying several times into the same ground state. The flatness of the distribution can be analytically derived by induction. (This can be shown by first noting that j N 0 0 i V decays with equal probability to j N

Four-level system
The situation considerably changes if we consider the full F g ¼ F e ¼ 1=2 four-level system, where je AE i V ≡ je AE1=2 i V . If we start again in je − i V , the single-particle decay rates are identical to Eq. (G2). However, a crucial difference from the three-level case is that now theΠ þΠ− andΣ þΣ− processes will lead to a temporary transfer of population into je þ i V .
In Fig. 17(c), we show that the ground-state distribution obtained in this case is different from the three-level case of Fig. 17(b) where je þ i V is not included. To analytically understand this, it is better to work in theẑ k basis where the initial state is ð1= ffiffi ffi 2 p Þðje − i k − je þ i k Þ [cf. Eq. (42)]. The superradiant decay withL − andR − essentially transforms the populations as je − i k → jg þ i k and je þ i k → jg − i k , but it kills all initial coherences between je − i V and je þ i V . Therefore, the final state isρ This corresponds to a ring of width ffiffiffiffi N p around the equator in a Bloch sphere where jg þ i k and jg − i k are the north and south poles defining the z axis. In this language, the probability distribution for the populations of jg AE i V is the projection of this ring onto the x axis. Interestingly, this leads to the same shape as the one obtained in Fig. 18 for the excited-state population when starting at a maximum in between two dark states.

APPENDIX H: DARK-STATE DISTRIBUTION OF EXCITATIONS
In this section, we analyze in closer detail the behavior of the six-level system of Sec. VIII for the case where the atoms start at a maximum of the superradiant potential VðθÞ. Specifically, we study the full distribution function of the excitation fraction n e that the system reaches as t → ∞.
In an experiment, this would correspond to the distribution of values of n e measured from shot to shot. To calculate this quantity, we simulate 10 5 TWA trajectories and compute the histogram of the values obtained that fall within bins of width 0.001. In Fig. 18, we show the results for a system starting at the θ 0 ¼ 2.822π maximum [cf. potential in Fig. 8(a)]. Interestingly, the distribution shows a doublepeak structure with a flat section in between. The peaks lie at n e ≈ 0.46 and n e ≈ 0.09, which correspond to the excited-state fractions of the dark states left and right from our initial θ 0 in the potential. As N is increased, we find that the peaks become sharper.
To test the validity of the TWA results in a dynamical regime driven by quantum fluctuations, we show in the inset of Fig. 18 ED results for up to N ¼ 30. The distribution obtained confirms the double-peak structure, although it is much smoother and slightly shifted compared to TWA. This behavior is, however, expected for small N. As N becomes larger, the peaks also get sharper and shifted in the ED simulations. While we cannot perform ED simulations for very large N, the tendency of the ED results is consistent with TWA.
The double-peak shape of the distribution in Fig. 18 resembles the histogram one would obtain from uniformly sampling a sin ϕ function within ϕ ∈ ½0; 2πÞ. The key to understanding this shape is that atoms starting at a maximum have zero mean dipole moment, D þ ¼ 0. In each TWA trajectory, however, the initial quantum fluctuations will lead to a small but nonzero dipole moment. For simplicity, this quantum noise can be parametrized by D þ ¼ ϵe iϕ with ϵ ≪ 1. This dipole essentially constitutes a small kick away from the maximum that leads to the decay of the atoms. Importantly, the direction of the dipole moment will be random as determined by ϕ. Because of this, in each TWA trajectory, the α-Bloch vectors will rotate in a different plane, which is equivalent to moving along a different potential landscape.
For a fluctuation with ϕ ¼ AEπ=2 (which implies noise in the S y α variables), the system will decay to the MF dark state on the right or left of θ 0 in Fig. 8(a). For any other ϕ, the system decays into a different dark state whose excitation fraction lies approximately in between the excitation values of the left (n D L e ) and right (n D R e ) dark states. Empirically and semianalytically, we then find that the final excited-state population of the MF dark state depends on ϕ approximately as n e ðϕÞ ¼ 1 2 ½ðn D R e − n D L e Þ sin ϕ þ ðn D R e þ n D L e Þ. For completeness, we mention that for states starting away from maxima the final distribution has a Gaussian shape, which is consistent with a MF coherent state.

APPENDIX I: FINDING TWO-POLARIZATION MEAN-FIELD DARK STATES
In general, we find that dark states are ubiquitous in the two-mode multilevel superradiance problem. When both polarizations are actively involved, making semianalytical predictions as for single-mode superradiance is hard. Even accurate dynamics simulations are nontrivial, as the discrepancies between the TWA and cumulant in Fig. 10 show. Nevertheless, it is still possible to numerically find twopolarization MF dark states in the following manner. At the mean-field level, dark states are simply states where the average dipole moment vanishes. If we make a product-state ansatz for the atoms as jDi ¼ jDi ⊗N i with jDi i ¼ P m α m jg m i q;i þ P n β n je n i q;i , we can find a MF dark state by requiring This condition translates into a small set of quadratic equations for α m and β n which can be solved numerically. Note that we can write Eq. (I1) with the two dipole moment operators corresponding to any polarization basis, e.g., hL AE i ¼ hR AE i ¼ 0. In order to assess whether the MF dark states found through Eq. (I1) are stable or unstable, we can then simply evaluate the second derivatives ðd 2 =dθ 2 ÞVðθÞ and ðd 2 =dθ 2 ÞUðθÞ for the desired state. Based on our numerical findings, we expect stable MF dark states to approximate quantum dark states as N → ∞.