Characterizing superradiant dynamics in atomic arrays via a cumulant expansion approach

Ordered atomic arrays with subwavelength lattice spacing emit light collectively. For fully inverted atomic arrays, this results in an initial burst of radiation and a fast build up of coherences between the atoms at initial times. Based on a cumulant expansion of the equations of motion, we derive exact analytical expressions for the emission properties and numerically analyze the full many-body problem resulting in the collective decay process for unprecedented system sizes of up to a few hundred atoms. We benchmark the cumulant expansion approach and show that it correctly captures the cooperative dynamics resulting in superradiance. For fully inverted arrays, this allows us to extract the scaling of the superradiant peak with particle number. For partially excited arrays where no coherences are shared among atoms, we also determine the critical number of excitations required for the emergence of superradiance in one- and two-dimensional geometries. In addition, we study the robustness of superradiance in the case of non-unit filling and position disorder.


I. INTRODUCTION
The interaction of dense atomic ensembles with light gives rise to a plethora of interesting many-body effects.One paradigmatic example is Dicke superradiance [1,2], a phenomenon in which the atoms in a totally inverted point-like sample synchronize and emit light coherently.This results in the cooperative speed-up of the atomic decay process and the emergence of a superradiant radiation burst [see Fig. 1(c)].Various aspects of Dicke superradiance have been studied [3][4][5][6][7][8] and experimentally observed for a broad variety of platforms, ranging from atomic gases [9][10][11][12][13][14][15] to solid state systems such as quantum dots [16,17], nitrogen-vacancy centres [18] and two-dimensional materials [19].
Over the last years, ordered atomic arrays with subwavelength lattice spacing have emerged as a promising platform to study collective light-matter coupling in free space.In the few excitation limit, these systems naturally exhibit collective superradiant and subradiant states [20].Various aspects of this low-excitation regime were studied in detail in recent years [20][21][22][23][24][25][26][27][28][29][30][31].However, analyzing the cooperative dynamics for large system sizes in the multi-excitation regime is challenging due to the rapidly growing Hilbert space.This restricts the study of ordered atomic ensembles with multiple excitations to very small systems of about ten atoms, for which the full solution of the open-system master equation is still feasible by means of Monte Carlo wave function (MCWF) methods [32][33][34].
For fully inverted atomic ensembles, it has been recently shown that the existence of a superradiant burst can be simply determined from the the statistics of the first two emitted photons [34][35][36][37].This method does not require to propagate the equations of motion in time, and therefore allows the efficient analysis of very large sys- (a) The EM vacuum mediates interactions among emitters trapped in a periodic array.Superradianct dynamics is characterized for fully or partially excited arrays where no coherences are shared among emitters at t = 0. (b) The complexity of the system can be reduced by truncating the number of equations of motion via a cumulant expansion of multi-order correlators.(c) This approach allows the characterization of superradiant dynamics for large system sizes.The bottom panel shows the total emission rate γtot/(N Γ0) of an examplary time evolution for a chain of N = 196 atoms and lattice spacings a = 0.3λ0 (solid line) and a = 0.5λ0 (dashed line).
tems (up to 10 6 atoms).It provides, however, no information about the decay dynamics, such as the magnitude of the superradiant peak or the existence of subradiance at late times.Alternatively, recent theoretical studies have developed an effective two-atom description of the many-body system capable of capturing superradiance in inverted three-and two-dimensional ensembles [38][39][40].
In this work, we employ cumulant expansions of the operator expectation values governing the system's dynamics to characterize the decay process of dipole-dipole coupled arrays of atoms [41][42][43][44].Based on neglecting high order quantum correlations, cumulant or cluster expansions drastically reduce the degrees of freedom used to describe the system.As opposed to MCWF methods [45,46], which require 2 N variables to describe an N -atom array, the cumulant expansion up to order n contains only ∼ N n terms (see sketch in Fig. 1).To this end, we derive the equations of motion for second and third order cumulant expansions in Section II, which provide analytical insights on the mechanisms leading to superradiance.We then benchmark the formalism in Section III by comparing the resulting dynamics with that predicted by MCWF.We demonstrate that cumulant expansions capture the early dynamics resulting in the superradiant burst with remarkable accuracy and find that they allow to simulate unprecedented system sizes of up to a few hundred atoms, more than one order of magnitude larger than with MCWF.
Based on these insights we then employ this formalism to unveil new insights on the physics of superradiant emission (section IV).In particular, we characterize the magnitude of the superradiant peak and analyze its scaling as a function of lattice spacing and atom number for large one-dimensional and two-dimensional inverted arrays in Section IV A. We additionally demonstrate that the enhanced atomic decay in inverted systems originates from a fast buildup of coherences among atoms.Hence, a natural question arises: what is the critical fraction of excited emitters required for this cooperative synchronization effect to occur if the atoms don't share coherences at initial times?In Section IV B, we analytically determine this critical condition for incoherently partially inverted arrays, and further compute the value of the superradiant peak versus excitation fraction using cumulant expansions.Finally, we extend the study to atomic arrays with finite filling fractions in Section IV C and with position disorder in Section IV D.

II. MODEL
We consider an ensemble of N identical two-level atoms that interacts with the vacuum electromagnetic field in free space.Applying the Markov approximation and integrating out the field degrees of freedom (see Fig. 1), one obtains the equations of motion for the atomic operators in the Heisenberg picture [47,48] Here, Ĥ and L( Ô) respectively describe the coherent and incoherent dipole-dipole coupling mediated by the vac- where σeg i = |e i g i | and σge i = |g i e i | are the raising and lowering operators for atom i, which is located at r i .Additionally, the coherent J ij and dissipative Γ ij couplings between atoms i and j are obtained from the Green's tensor for a point dipole in vacuum G, given in Appendix A, as where ω 0 = 2πc/λ 0 corresponds to the transition wavelength of the emitters, d to their transition dipole moment and r ij = r i − r j to the vector connecting atoms i and j.Here, Γ ij = Γ 0 is given by the spontaneous decay rate of a single atom in vacuum and the Lamb-shift J ii is included in the definition of ω 0 .For the remainder of this work, we consider a transition dipole moment perpendicular to the array.Note that similar results are obtained for other polarizations.
Finally, the last term in Eq. (1), F( Ô), describes the quantum Langevin noise arising from vacuum fluctuations [47].Assuming white noise for F( Ô), the expectation value F( Ô) vanishes.Since we are ultimately interested in averages over atomic operators Ô , we drop this term in the following discussion to simplify notation.
In this work, we consider initial states where N exc atoms are incoherently excited at t = 0, where |g is the state with all atoms in the ground state and E denotes the set of initially excited atoms.The two particle correlations σeg i σge j vanish for these states.Hence, no coherences are shared among atoms at initial times even for partially excited arrays.Such states can be prepared by either using a spatially incoherent light source to excite the atomic atomic array or by imposing a random detuning pattern over the duration of the excitation pulse.
The equations of motion of the relevant first-and second-order operators σee i , σeg i σge j and σee i σee j can be readily obtained from Eq. ( 1) Operators of order n, i. e., involving n products of individual atomic operators, are generally coupled to operators of order n + 1.This leads to a coupled set of differential equations for system operators up to order N , where N denotes the particle number.

A. Cumulant expansion
The number of equations needed to fully describe the system grows exponentially with atom number, limiting the system sizes that can be numerically simulated to about sixteen atoms.To study larger systems, one can take expectation values of the equations of motion in Eq. (6) and truncate the set of equations by approximating averages of higher-order operators with combinations of averages of lower-order ones (see Fig. 1).This method is known as cumulant or cluster expansion [41][42][43][44].For initial states given by Eq. ( 5), only the three second-order operators given in Eq. ( 6) develop non-zero expectation values σee i , σee A similar procedure can be used to derive the thirdorder cumulant expansion of the atomic system, which additionally includes the expectation values σee

B. Emission properties
While a recent theoretical work [44] focused on employing a cumulant expansion approach to determine twotime correlation functions, we hereby focus on characterizing the emission properties of the atomic array.In fact, the equations of motion corresponding to the sec- The inset in panel (c) represents the total emission rate at late times in logarithmic scale.The atoms are considered to be polarized in the direction perpendicular to the lattice plane throughout this work.
ond order cumulant expansion provide a powerful tool to elucidate the fundamental mechanisms resulting in superradiance.
The decay process of the atomic ensemble can be characterized by the excited-state population and the total emission rate of the system Eq. ( 8) provides a very nice intuitive picture on the underlying mechanisms resulting in superradiance.For independent particles (Γ ij = 0) the emission rate is proportional to the individual atom decay rate Γ 0 .For interacting particles, however, the term ∝ Γ ij σeg i σge j modifies the emission rate, and can eventually result in a superradiant burst if the interactions Γ ij and pair correlations σeg i σge j are sufficiently large.The existence of a radiation burst can be predicted from the derivative of the total emission rate at t = 0 [37].For initial states given by Eq. ( 5), we obtain where Ô 0 ≡ Ô (t = 0).For γtot,0 > 0, the emission rate initially increases and the dynamics consequently result in a superradiant peak.For γtot,0 < 0, on the other hand, γ tot initially decreases, indicating that the dissipative interactions between emitters cannot build the amount of coherences required for the appearance of a burst.For an initially fully inverted system with σee i 0 = σee i σee j 0 = 1, Eq. ( 9) reduces to and is identical to the expressions derived in Ref. [37], as well as to the expression obtained via the two-photon correlation function in Ref. [34].Note that both Eqs. ( 9) and ( 10) are exact, as second-order cumulants perfectly capture two-photon processes at initial times.While typically sufficient when studying the decay of fully inverted systems or partially excited coherent spinstates [34], this criterion can fail to identify radiation bursts for more general initial conditions.One example are incoherently driven, partially excited arrays, whose initial states are in a superposition of dark and bright states as it is presented in Ref. [49].In this case, strong coherent dipole-dipole interactions J ij can generate a net population transfer from subradiant to superradiant states, which ultimately leads to This phenomenon, which is triggered by the Hamiltonian evolution of the system, only shows signatures in the emission rate for t > 0. It can consequently not be captured by Eq. ( 9), which solely relies on the dissipative couplings and the atomic properties at zero time.

III. BENCHMARKING CUMULANT EXPANSIONS
When simulating the dynamics of coupled spin systems, there is always a trade-off between the maximum number of particles that can be numerically simulated and the accuracy of the resulting dynamics.Increasing the order of the cumulant expansion, for example, reduces the former and generally improves the latter.In this section, we study which expansion order is required to correctly capture the evolution of the atomic ensemble.This analysis can only be performed for small system sizes, for which a full solution of the master equation is possible.
In Fig. 2, we plot dynamic properties of a fully inverted ten-atom chain with lattice spacing a = 0.1λ 0 using first (grey dotted line), second (blue dash-dotted line) and third (red solid line) order cumulant expansions, as well as the solution of the full master equation (black dashed).Third order cumulants exhibit very good agreement with the exact time evolution and properly capture the magnitude of the superradiant peak (maximum of γ tot /N ), while second order cumulants slightly overestimate it.The outburst of radiation is accompanied by a buildup of the atomic coherences σeg i σge j at early times due to the repeated application of the same bright jump operators [34].At late times, however, the performance of the cumulant expansion worsens.In particular, the second order cumulant expansion fails to correclty predict the late dynamics [see inset in Fig. 2(c)], and the third order cumulant expansion sometimes results in unphysical behaviors such as growing excitation number in the absence of drive [44].This poor performance at late times can be attributed to the population of subradiant states during the decay process.These states are typically highly entangled, and therefore contain many high-order inter-atom coherences that lead to destructive interference of the electromagnetic field emitted by different atoms.Consequently, an expansion that neglects higher-order coherences performs worse in predicting these quantities.For very large system sizes (N > 100) and small lattice spacings, the cumulant expansion can also result in unphysical behaviors such as a growing excited state population in the absence of drive.
Finally, it is worth noting that the first-order cumulant expansion, commonly referred to as the mean-field approximation, simply results in an independent exponential decay of the atomic ensemble (without a buildup of the atomic coherences) and does not capture neither super-nor subradiance.
To benchmark the accuracy of cumulant expansions of different orders, we compare two figures of merit of the decay of inverted atomic arrays as a function of particle number N and lattice spacing a.In Fig. 3(a), we plot the maximum emission rate or magnitude of the superradiant peak, max(γ tot /N ), as a function of particle number N for a one-dimensional chain with lattice spacing a = 0.1λ 0 .The inset shows the error made by second-order (blue circles) and third-order (red diamonds) cumulant expansions, defined as ∆ peak master = max(γ cum tot )/ max(γ master tot ) − 1, in percent.Not surprisingly, third-order cumulant expansions are more accurate than second-order ones and, more importantly, their error does not grow linearly with N .Fig. 3(b) demonstrates that the accuracy considerably worsens with decreasing lattice spacing.This effect can be attributed to the rapidly growing energy shifts, which generate correlations between atoms.
Another relevant quantity that characterizes the decay process of the ensemble is the subradiant population p sub , that is, the amount of excitation left in the atomic array once the system enters the subradiant regime.This occurs once the instantaneous decay rate γ inst ≡ γ tot /(p exc Γ 0 ) goes below 0.1, as defined in Ref. [49].Again, third-order cumulants provide an excellent estimate of this quantity, while second-order cumulants fail to capture it.The larger errors of secondorder cumulants when estimating p sub as opposed to max(γ tot /N ) are due to the fact that the former contains information of the decay process until times much longer than the radiation burst, where subradiant evolution starts to dominate.Still, p sub largely depends on the superradiant dynamics and is therefore accurately captured by the third-order cumulant expansions [49], even if the subsequent subradiant evolution is not [see inset of Fig. 2(c)].

IV. SUPERRADIANCE
Having elucidated the accuracy with which second and third order cumulant expansions capture superradiance, we now investigate the conditions under which this phenomenon emerges in one-and two-dimensional atomic arrays, as well as the magnitude of the radiation burst reached by the ensemble of emitters.

A. Fully inverted arrays
Superradiant emission is typically studied in fully inverted arrays, where all atoms have been initially excited by an intense laser drive.In this scenario, the existence of a burst as a function of system size and lattice spacing can be inferred from Eq. ( 10) [34][35][36][37].This condition is only based on the structure of the equations of motion at t = 0 and no predictions on the resultant time evolution can be made.Here, we use the cumulant expansions approach outlined above to extend previous studies [34-37, 39, 40, 44] on superradiance in ordered subwavelength arrays by analyzing the time evolution leading to superradiance.This allows us to quantify the magnitude of the superradiant peak for system sizes of a few hundred atoms, compatible with state of the art experiments.
In Fig. 4 we determine the superradiant regime for different particle numbers and lattice spacings for a chain of atoms [Fig.4(a)] and a square lattice [Fig.4(b)].We calculate the time evolution for a certain set of parameters (N and a) by solving the cumulant equations up to second order and determine the maximum of the total emission rate, which is plotted as a color code in Fig. 4(a) and (b).A superradiant peak occurs if max(γ tot ) > N γ 0 , that is, if the maximum emission rate is larger than that of independent decay.The white solid line separates the regime where a superradiant burst occurs from the regime where no cooperative enhancement of the atomic emission takes place, and coincides with the critical value obtained from Eq. (10).The results presented in Fig. 4(a) and (b) demonstrate that the transition between both regimes is not sharp and suggest that lattice spacings well below the ones predicted in Ref. [34][35][36][37] [white line Fig. 4(a) and (b)] are required to experimentally observe the initial speed-up of radiation.
Fig. 4(c) and (d) show the scaling of the superradiant peak with particle number for a 1D chain and a square lattice, respectively ( obtained with a third order cumulant expansion).We find that the maximum emission rate saturates for N → ∞ for a one-dimensional chain, whereas it increases as ∝ N β for a two-dimensional configuration.Performing linear fits of the curves shown in Fig. 4(d), we find that the exponent β increases for decreasing lattice spacing until it saturates for a/λ 0 → 0 [see Fig. 4(e)].It is worth noting that the values of β found in this work are consistent with those reported in Ref. [39].

B. Partially inverted arrays
For fully inverted arrays, the mechanism resulting in superradiance relies on the fast build up of coherences between pairs of atoms [see Eq. ( 8) and Fig. 2(b)].The natural question arises: how many atoms need to be excited for these coherences to build up fast enough and for a superradiant burst to occur?For that, we consider atomic ensembles driven by an incoherent driving field, such that only N exc randomly selected emitters are excited.This type of initial state, where σeg i σge j 0 = 0, can be generally described by Eq. ( 5).In the inset of Fig. 5(a), we show in grey the emission rate for different initialization of a square lattice with thirty-six atoms, thirty of which are initially excited.While the exact magnitude of the burst depends on the specific configuration of excited atoms, one can extract a characteristic peak size by averaging over all trajectories (thick blue line).The resulting average peak size is then plotted as a function of the fraction of excited atoms, n exc := N exc /N , in Fig. 5(a).We find that there is a critical excitation fraction, n crit exc , required to obtain a superradiant burst both for atomic chains and two-dimensional arrays, which depends both on system size and the specific geometry.
One can analytically compute n crit exc by averaging Eq. ( 9) over all possible configurations of excited atoms.Defining N de = N − N exc as the number of de-excited atoms, we obtain an average derivative of the emission rate at t = 0 equal to (for details see Appendix D) A superradiant burst is observed for γtot,0 > 0, that is, for This analytical result shows perfect agreement with the numerically obtained results [see vertical dashed and dash-dotted lines in Fig. 5(a)].
In the Dicke limit a → 0, where all emitters are located at the same spatial position, all dissipative couplings are equal to the spontaneous decay rate, i. e., Γ mn = Γ 0 ∀m, n, and Eq. ( 12) results in n crit exc = 1/2 + 1/N .That is, an excitation fraction larger than one half, i. e., N exc > N/2 + 1, is required to observe a burst.This can be intuitively understood by realizing that the symmetric Dicke state with maximum decay rate corresponds to the state where half of the atoms are excited.Hence, an initially fully inverted system dynamically evolves into states with higher and higher decay rates while cascading down the Dicke ladder, which ultimately results in the appearance of a superradiant emission peak.Once half of the atomis are de-excited this trend reverses and the system starts populating states with decreasing decay rates, causing the peak to vanish.If only half (or less) of the atoms are excited initially, the system will never decay into states with a larger decay rate than the initial one, and the atomic ensemble is bound to decay exponentially without the emergence of a burst.
The analytical result in Eq. ( 12) provides insight into the role of the lattice geometry and dimensionality for arbitrary system sizes.For extended arrays of arbitrary dimension, the sum in the denominator of Eq. ( 12) is always smaller than for the Dicke case, and a larger excitation fraction is needed to attain superradiance.In one-dimensional arrays, the critical excitation fraction tends to a constant for large particle number, as shown in Fig. 5(b).In two-dimensional arrays [see Fig. 5(c)], however, it decays logarithmically with system size N , lim N →∞ n crit exc /N ∼ 1/2 + A/ ln( √ N ), and eventually reaches the Dicke limit for infinite systems.For completeness, it can be shown that this scaling is improved in three-dimensional lattices, where lim N →∞ n crit exc /N ∼ 1/2 + A × N −1/3 [36,37].Figs.5(b) and (c) also show that a larger excitation fraction is needed to attain superradiance in atomic chains than in a square lattice for equal atom number N and lattice spacing a, as well as the fact that smaller values of a result in a smaller n crit exc .Hence, Fig. 5 and Eq.(12) show that there exists a well-defined critical light intensity of the excitation pulse Magenta diamonds correspond to a one dimensional chain and cyan diamonds to a two dimensional square lattice, both with spacing a/λ0 = 0.1.The inset shows the time evolution of 100 trajectories, each of them corresponding to a different, random initial distribution of thirty excitations.The blue curve indicates the average over all trajectories and is used to determine the magnitude of the superradiant peak shown in the main figure.The vertical dashed and dashdotted lines correspond to the critical excitation fraction, determined by Eq. ( 12) for atomic chains and two-dimensional lattices, respectively.(b)-(c) Critical excitation, obtained via Eq.( 12), as a function of particle number for (b) onedimensional and (c) two-dimensional arrays with different lattice spacings.For a/λ0 → 0 the system approaches the Dicke case (black dashed line).The blue circles (a/λ0 = 0.1), red diamonds (a/λ0 = 0.15), yellow crosses (a/λ0 = 0.2) and green hexagons (a/λ0 = 0.25) are obtained numerically and show a good agreement with the analytical result.
for which superradiance emerges.It is worth noting, however, that Eq. ( 12) can fail to identify superradiant bursts triggered by the Hamiltonian dynamics at times later than t = 0 (for details see Refs.[38,49]).While these may result in superradiant peaks appearing for n exc < 0.5, they require strong coherent couplings J ij or large dephasing or Doppler broadening and therefore do not typically appear for lattice spacings a ≥ 0.1λ 0 considered in Fig. 5.
Finally, it is worth noting that, in experimental implementations of superradiance, the excitation pulse might not be able to fully invert the atom system.Our results show that superradiance is reasonably robust to finite excitation fraction, as a radiation outburst still occurs at low enough lattice spacings even if only 90% or less of the atoms are initially excited.14), as a function of particle number for (b) one-dimensional and (c) two-dimensional arrays with different lattice spacings.For a/λ0 → 0, the system approaches the Dicke case (black dashed line).The blue circles (a/λ0 = 0.1), red diamonds (a/λ0 = 0.15), yellow crosses (a/λ0 = 0.2) and green hexagons (a/λ0 = 0.25) are obtained numerically and show a good agreement with the analytical result.

C. Finite filling fractions
Typically, experimental implementations of atomic arrays also suffer from a finite filling fraction of the lattice.In Ref. [21], for example, which demonstrated that twodimensional array with subwavelength spacing can reflect incident beams of low intensities, approximate filling fractions of ninety percent were achieved.Here, we show that superradiance for fully inverted arrays is very robust to missing atoms and that large superradiant bursts can still be achieved for filling fractions much lower than those attainable in state of the art experiments.
Fig. 6(a) shows a similar analysis as it was performed in the previous section.Now, however, all atoms are initially excited and only a fraction η ≡ N filled /N of the N available lattice sites is populated by an excited atom.In other words, the lattice contains N hol = N − N filled holes or missing atoms.Fig. 6 shows the average magnitude of the peak for one-and two-dimensional arrays with N = 36 lattice sites as a function of the filling fraction η.The peak values are again obtained by averaging over time evolutions for 100 random distributions of holes.The magnitude of the peak slowly decreases with diminishing η, and finally reaches zero at a critical filling fraction ñcrit exc that depends on the properties of the atomic ensemble.Again, we can estimate η crit by averaging γ0 over all possible configurations of missing atoms (for details see On average, a burst occurs for γtot,0 > 0, that is, for In the Dicke limit, Eq. ( 14) reduces to the well known expression η > 2/N , i. e., N filled > 2. For two atoms, the decay rates from |ee to (|eg + |ge )/ √ 2 and from (|eg + |ge )/ √ 2 to |gg are both 2Γ 0 .The two-atom system therefore decays faster than in vacuum, but cannot develop a superradiant burst.For this to occur, the system needs to dynamically access states with larger and larger decay rates while cascading down the Dicke ladder, which naturally occurs for ensembles with three or more atoms.Additionally, it is straightforward to check that n crit exc = 1/2 + η crit /2.Thus, the scalings presented in Section IV B also hold for the case of missing atoms.
It is worth noting that the excitation fraction needed to observe a peak is considerably larger in the case of partially excited configurations, presented in Section IV B, than in the case of missing atoms or holes.The main difference between both scenarios consists on the fact that the latter remains a fully inverted system, while the former does not.That is, an inverted array with missing atoms can simply cascade down the whole ladder of superradiant states by repeated action of the brightest jump operator -which is now weaker than in a perfect array-and can therefore attain a burst even for very small filling fractions.However, this is not the case for incoherent partially inverted arrays, whose initial state is already a superposition of bright and dark states in an intermediate sector of the state space (see also Ref. [49]) and therefore has a reduced likelihood to quickly build up the necessary correlations to create a superradiant burst.

D. Position disorder
Another potential imperfection that might influence the superradiant burst in realistic setups is disorder in the atomic positions.This can either be due to some distortion in the underlying trapping potential or due to atomic motion in the respective trapping potentials.
To take position disorder into account, we sample the modified atomic positions from a normal distribution N d (r 0 , σ), where σ denotes its standard deviation and r 0 the positions of the disorder-less array.We then average over 100 different excitation and disorder patterns and plot the determined peak height in Fig. 7.We find that the effect of disorder remains small as long as σ is below 0.1a/λ 0 , but rapidly increases for larger σ (see green curves in Fig. 7).Still, these results suggest that superradiance is very robust against position disorder.

V. CONCLUSIONS
We performed a numerical and analytical in-depth analysis of superradiance in ordered atomic arrays based on a cumulant expansion approach.This formalism allowed us to simulate the dynamics of open quantum systems for large particle numbers and to gain insights on the physics behind super-and subradiance.In particular, we identify the scaling of the superradiant peak with particle number both for one-dimensional and twodimensional arrays.Additionally, we show that there exists a critical excitation fraction above which a superradiant burst occurs, and demonstrate that superradiance is a robust phenomenon that prevails in the presence of position disorder and finite filling and excitation fraction.
The results presented here are expected to be readily observable in various state-of-the-art platforms, ranging from atomic tweezers [50] and optical lattices of cold atoms [21] to solid state platforms such as twodimensional materials [19,51] or vacancy centers in crystals [18].Our work also paves the way towards efficiently calculating dynamic light emission patterns from atomic arrays, as well as towards devising novel schemes to prepare multi-excitation subradiant states [49].Finally, the existence of a critical excitation fraction to attain superradiance exhibits a certain resemblance with the superradiant transition in the driven dissipative Dicke model [52][53][54][55].Exploring this connection is an exiting avenue for future work.The Green's function for a point dipole in free space used in Eq. ( 4) can be written in Cartesian coordinates as [56,57] where k = ω/c, r = |r|, and α, β = x, y, z.
Appendix B: Second-order cumulant expansion For initial states with no correlations as defined in Eq. ( 5), the only nonzero expectation values of first-and second-order operators at t = 0 are σee Using Eq. ( 1), one can show that the only additional expectation values that become non-zero during the evolution of the system are σeg i σge j .This confirms that Eq. ( 6) is sufficient to describe the dynamics of the atomic ensemble up to second order.
Taking expectation values in Eq. ( 6) and replacing averages over third-order operators by [41,43] Appendix C: Third-order cumulant expansion The third-order cumulant expansion is obtained by deriving the Heisenberg equations of motion for operators up to third order and replacing the averages of fourth-order operators by lower-order ones via the rule [41,43] For the initial states given in Eq. ( 5), the non-zero expectation values at t = 0 are those given in Eq.( B1) plus where (a ↔ b) indicates that an additional term appears equal to the previous one but with indexes a and b swapped.

Appendix D: Critical fraction of excited atoms
Let us consider an ensemble of N atoms.At initial timex, N exc of them are excited, while the remaining N de = N − N exc are de-excited (see Section IV B).We label as E and D the set of excited and de-excited atoms, respectively.We consider incoherently excited arrays, such that σee i 0 = 1 if i ∈ E and zero otherwise, and σee i σj 0 = 1 if i, j ∈ E and zero otherwise.The derivative of the total emission rate for any configuration, given in Eq. ( 9), can be expressed as γtot,0 = − We further average over all possible configurations of N de de-excited atoms.We label the de-excited atoms with indexes α 1 to α N de , each of them taking values from 1 to N .Noting that no two indexes can be equal, the total number of permutations is N !/(N − N de )!.Then, the average derivative of the total decay rate can be written as γtot,0 = −Γ 2 0 N exc + i j =i where Σ 1 and Σ 2 are        Note that the population at the lattice site of a missing atom is always zero, and so is its second derivative too.
The population in a filled lattice site, on the other hand, is equal to that of an alternative, fully inverted system containing only N filled atoms at the occupied positions.Then, γtot,0 can be written as γtot,0 = −Γ 2 0 N filled + i∈F j∈F ,j =i which takes the same form as Eq.(D5) but with different pre-factors.Using Eq. (D7) and Eq.(D8), one readily finds the average derivative of the total emission rate, given in Eq. ( 13) of the main text.

Figure 1 .
Figure 1.Schematic illustration of the theoretical approach.(a)The EM vacuum mediates interactions among emitters trapped in a periodic array.Superradianct dynamics is characterized for fully or partially excited arrays where no coherences are shared among emitters at t = 0. (b) The complexity of the system can be reduced by truncating the number of equations of motion via a cumulant expansion of multi-order correlators.(c) This approach allows the characterization of superradiant dynamics for large system sizes.The bottom panel shows the total emission rate γtot/(N Γ0) of an examplary time evolution for a chain of N = 196 atoms and lattice spacings a = 0.3λ0 (solid line) and a = 0.5λ0 (dashed line).
evolution.Approximating the three-operator averages in Eq. (6) as σee a obtains a closed set of differential equations for the first-and second-order expectation values, which correspond to the second-order cumulant expansion (see Appendix B).
explicit form of the corresponding equations of motion is presented in Appendix C.

Figure 2 .
Figure 2. Comparison between the dynamics obtained via cumulant expansion of orders 1-3 and the solution of the full master equation for a chain of N = 10 atoms with lattice spacing a = 0.1λ0.(a) Excited state population pexc(t), (b) average pair correlations and (c) total emission rate γtot/(N Γ0).The inset in panel (c) represents the total emission rate at late times in logarithmic scale.The atoms are considered to be polarized in the direction perpendicular to the lattice plane throughout this work.

Figure 3 .
Figure 3. Benchmarking the cumulant expansion method.(a) Maximum value of the superradiant peak for a fully inverted chain of atoms with lattice spacing a = 0.1λ0 as a function of atom number.(b) Value of the superradiant peak for a N = 10 atom chain as a function of lattice spacing.The insets in panels (a) and (b) show the relative error between the exact master equation solution and the cumulant expansions in percent.(c)-(d) Magnitude of the subradiant population p sub , i. e., the excited population left in the array by the time the instantaneous decay rate γinst = γtot/(pexcΓ0) = 0.1, as a function of (c) particle number and (d) lattice spacing.Same parameters as in panels (a) and (b).

Figure 4 .
Figure 4. (a)-(b) Maximum emission rate max(γtot/N Γ0) for (a) an atomic chain and (b) a two-dimensional square lattice as a function of particle number N and spacing a/λ0 obtained via a second order cumulant expansion.The white solid line separates the region where a superradiant burst is observed from the region where no superradiance occurs.(c)-(d) Scaling of the superradiant peak as a function of particle number for (c) a chain and (d) a square lattice for different lattice spacings.(e) Exponent β characterizing the power law dependence of the superradiant peak with particle number, ∝ N β , for a square lattice.The values are obtained from the linear fits in panel (d), shown as dash dotted lines.

Figure 5 .
Figure 5. (a) Value of the superradiant peak as a function of the initial incoherent excitation fraction for N = 36 atoms.Magenta diamonds correspond to a one dimensional chain and cyan diamonds to a two dimensional square lattice, both with spacing a/λ0 = 0.1.The inset shows the time evolution of 100 trajectories, each of them corresponding to a different, random initial distribution of thirty excitations.The blue curve indicates the average over all trajectories and is used to determine the magnitude of the superradiant peak shown in the main figure.The vertical dashed and dashdotted lines correspond to the critical excitation fraction, determined by Eq. (12) for atomic chains and two-dimensional lattices, respectively.(b)-(c) Critical excitation, obtained via Eq.(12), as a function of particle number for (b) onedimensional and (c) two-dimensional arrays with different lattice spacings.For a/λ0 → 0 the system approaches the Dicke case (black dashed line).The blue circles (a/λ0 = 0.1), red diamonds (a/λ0 = 0.15), yellow crosses (a/λ0 = 0.2) and green hexagons (a/λ0 = 0.25) are obtained numerically and show a good agreement with the analytical result.

Figure 6 .
Figure 6.(a) Value of the superradiant peak for finite filling fractions η for a fully inverted array.Magenta diamonds correspond to a one dimensional chain and cyan diamonds to a two dimensional square lattice, both with spacing a/λ0 = 0.1.(b)-(c) Critical filling fraction to observe a superradiant peak, obtained via Eq.(14), as a function of particle number for (b) one-dimensional and (c) two-dimensional arrays with different lattice spacings.For a/λ0 → 0, the system approaches the Dicke case (black dashed line).The blue circles (a/λ0 = 0.1), red diamonds (a/λ0 = 0.15), yellow crosses (a/λ0 = 0.2) and green hexagons (a/λ0 = 0.25) are obtained numerically and show a good agreement with the analytical result.

Figure 7 .
Figure 7. Effect of position disorder, sampled from a Gaussian distribution with standard deviation σ = (0.01, 0.05, 0.1, 0.2)a (blue circles, red diamonds, orange crosses, green squares) for a chain (a) and a square lattice (b) of N = 16 atoms and averaged over 50 different lattice configurations and initial excitation distributions.The dashed lines show the result for a perfect lattice without position disorder.The insets exhibit the difference ∆ = max{γtot/(N Γ0)}|σ=0 − max{γtot/(N Γ0)}|σ between the zero disorder case and the disordered systems.Disorder affects superradiance if σ > 0.1a.
atom i and j are excited 0, otherwise.

one obtains a closed set of differential equations for
atom i, j and k are excited 0, otherwise.(C2)Again,one can show that the the decay process only couples these expectation values of populations with σeg the expansion in Eq. (C1), one finally obtains the third-order cumulant expansion
and finally result in Eq. (11) of the main text.
Appendix E: Critical filling fractionLet us consider an array with N lattice position, N filled of which are occupied and N hol = N − N filled are empty (see Section IV C).Labelling F and H as the set of filled and empty lattice sites, we can compute the derivative of the total emission rate in a similar manner as in Appendix D. Now, the second order derivatives of the population in lattice site i read

ACKNOWLEDGMENTS
We would like to thank Valentin Walther and Yidan Wang for fruitful discussions.O.R.B. acknowledges support from Fundación Mauricio y Carlota Botton and from Fundació Bancaria "la Caixa" (LCF/BQ/AA18/11680093).S.O. is supported by a postdoctoral fellowship of the Max Planck Harvard Research Center for Quantum Optics.SFY would like to acknowledge funding from NSF through the CUA PFC and the QSense QLCI as well as from AFOSR.O.R.B. and S.O.contributed equally to this work.