Breakdown of steady-state superradiance in extended driven atomic arrays

Recent advances in generating well controlled dense arrangements of individual atoms in free space have generated interest in understanding how the extended nature of these systems influences superradiance phenomena. Here, we provide an in-depth analysis on how space-dependent light-shifts and decay rates induced by dipole-dipole interactions modify the steady-state properties of coherently driven arrays of quantum emitters. We characterize the steady-state phase diagram, with particular focus on the radiative properties in the steady-state. Interestingly, we find that diverging from the well-established Dicke paradigm of equal all-to-all interactions significantly modifies the emission properties. In particular, the prominent quadratic scaling of the radiated light intensity with particle number in the steady-state -- a hallmark of steady-state Dicke superradiance -- is entirely suppressed, resulting in only linear scaling with particle number. We show that this breakdown of steady-state superradiance occurs due to the emergence of additional dissipation channels that populate not only superradiant states but also subradiant ones. The additional contribution of subradiant dark states in the dynamics leads to a divergence in the time scales needed to achieve steady-states. Building on this, we further show that measurements taken at finite times for extended atom ensembles reveal properties closely mirroring the idealized Dicke scenario.


I. INTRODUCTION
Dicke superradiance, i. e., the effect where a group of excited atoms collectively emit radiation faster than they would individually, is a fundamental effect in light-matter interacting systems [1].Since Dicke introduced the concept of cooperative radiation of atomic ensembles [2], a series of seminal works have investigated cooperative phenomena in systems with strong light-matter interaction in various settings.Prominent examples include Dicke superradiance of atomic ensembles coupled to optical resonators and waveguides, in large clouds of Rydberg atoms and dense ensembles of two-level quantum emitters [3][4][5][6][7][8][9].More recently the advent of novel technologies which offer more control over the arrangement of individual atoms [10][11][12][13][14][15][16] or quantum emitters [17][18][19][20] has motivated an in-depth study of the role of the emitter geometry on superradiant emission [21][22][23][24][25][26][27][28][29].While the dissipative long-range nature of the underlying spin model makes exact studies of these systems challenging, it also provides opportunities to unveil fundamental governing principles of dissipative many-body systems, and to develop applications for novel lasing technologies [30], the preparation of entangled states [31,32] or the generation of non-classical states of light [33].
Most works characterizing superradiance focus on the transient dynamics of fully inverted ensembles.Recently, however, a seminal experimental work investigated the steady-state properties of a dissipative cooperative system under constant classical laser drive [34].In the Dicke limit, where all atoms are assumed to be located at the same spatial position, and the dissipative interactions among them consequently have the same strength, this * stefanostermann@g.harvard.eduFigure 1.Sketch of the considered model.A dense periodic array of two-level quantum emitters is driven by an external laser with Rabi frequency Ω.If the lattice spacing a is much smaller than the transition wavelength λ0 = 2πc/ω0, lightinduced dipole-dipole interactions will result in non-trivial collective dissipative dynamics.
system is known to exhibit quadratic scaling of the emitted light intensity with particle number in the steadystate [35][36][37] for sufficiently strong laser drive.For extended ensembles of atoms, geometry dependent energy shifts and decay rates are expected to alter these dynamics.
In this work, we investigate the role of light-induced cooperative shifts and decay rates for steady-state superradiance in driven-dissipative periodic arrays of quantum emitters under strong drive.While previous work focused on mapping out the steady-state phase diagram of this system [33], or perturbative approaches in the weak driving regime for disordered clouds [38,39], our focus lies on the full characterization of the superradiant emission properties as a function of the system parameters.The core finding of this work is that going beyond the paradigmatic Dicke limit modifies the radiative steadystate properties significantly.While the radiated light intensity in the steady-state scales quadratically with system size in the Dicke limit, it only scales linearly for extended ensembles (see Sec. IV).For very dense extended ensembles measured at finite times, however, we show that the scaling approaches that of the Dicke limit again, namely N 2 (see Sec. V).

II. MODEL
We model the system depicted in Fig. 1(a) as a set of two-level atoms with ground state |g i ⟩ and excited state |e i ⟩, located at positions r i in the x-y plane with a transition dipole matrix element d.To reduce the number of tunable parameters throughout this work, we choose d = (0, 0, 1) T .Analogous results, however, can be found for other polarizations.The Hamiltonian is then given as (we set ℏ = 1 here and for the remainder of this work) Ĥ = Ĥint + Ĥdrive , with where is the raising (lowering) operator for atom i, ∆ = ω − ω 0 denotes the detuning between the drive laser frequency ω and the atomic transition frequency ω 0 (including the Lamb shift J ii ), and J ij is the coherent interaction strength between distant emitters.For the analysis below we choose ∆ = 0 unless stated otherwise.The ensemble is driven with a plane wave drive with constant Rabi frequency Ω and wavevector k.The direction of k, i. e., the direction of the incoming pump beam, does not alter the steady-state properties, but it affects the transient dynamics and therefore is an important quantity to consider when analyzing finite time measurements in Sec.V.The strength of the coherent interaction is determined via [40,41], where G(r, ω) is the Green's tensor for a point dipole in vacuum [42,43] given in Appendix A, and r ij = r i − r j is the vector connecting atoms i and j.
The dissipative nature of the system is described by a Lindbladian of the form where The full system dynamics is then governed by the master equation for the atomic density matrix ρ [40,41,44]  Notably, the resulting long-range interacting spin model is not integrable, and the size of the density matrix ρ scales exponentially in system size.While a truncation of the Hilbert space into the so-called single excitation subspace is feasible in the weak driving regime, studying the effects of a strong drive requires to take into account the full exponential Hilbert space.This makes an exact quantum mechanical treatment of this many-body problem impossible for large system sizes, and has traditionally rendered the study of dipole-coupled ensembles in the strong driving regime unfeasible.To circumvent this limitation, we use approximative tools (mean-field approximation and second order cumulant expansions [45][46][47][48] in the analysis below to model larger systems and justify our findings by providing intuition based on the full quantum model in Eq. (3) for small atom numbers.
An important limit of this model is the so-called Dicke limit, which corresponds to N indistiguishable atoms lo-cated at a single position, for example r i = 0.The governing equations can be obtained from Eqs. (1a) and ( 2) by setting J ij = 0 and Γ ij = γ 0 , as well as defining collective raising and lowering operators Ŝ± = N i=1 σ± i .If the system is initialized in the ground state, the collective dynamics is restricted to the fully symmetric spin sector |S = N/2, m = −N/2...N/2⟩, often also referred to as the Dicke ladder.The Hamiltonian for the coherent drive then reduces to H Dicke drive = Ω/2( Ŝ+ + Ŝ− ) and the Lindbladian can be written as L Dicke The fact that the restriction to the symmetric subsector no longer holds in the free space case for finite lattice spacings is pivotal for the results discussed below.

III. STEADY-STATE PHASE DIAGRAM
We first determine the steady-state phase diagram as a function of Ω and a.It is convenient to express the Hamiltonian Ĥ in the Pauli basis consisting of the Pauli matrices σx,y,z by employing the relation σ± i = (σ i x ± iσ i y )/2.Then, the definition of the collective spin operators Ŝx,y,z = 1/2 N i=1 σx,y,z allows the interpretation of the results based on a collective spin on a single Bloch sphere.An exemplary phase diagram for a square consisting of four atoms is shown in Fig. 2(a).It exhibits two distinct regions indicated in blue and red.The blue region corresponds to the trivial magnetized phase where all the spins are polarized, whereas the red region marks the paramagnetic phase where ⟨ Ŝz ⟩ ≈ 0. Note that related characterizations of the superradiant phase transition in the Dicke limit were presented in Refs.[34][35][36][37].In the Dicke limit, the saturated paramagnetic regime also corresponds to the regime exhibiting enhanced superradiant emission, as will be further discussed in Sec.IV.In the Dicke case, the transition between the magnetized and the paramagnetic phase can be inferred by a peak in the value for the y component of the collective spin ⟨ Ŝy ⟩ [see Fig. 3(a)].In the finite size free space case, this peak is no longer as pronounced [see Fig. 3(b)].Nevertheless, the maximum of ⟨ Ŝy ⟩ remains a good indicator for the transition point between the magnetized and paramagnetic phase.Hence, we use max{⟨ Ŝy ⟩} to characterize the transition between these two regimes for the remainder of this work.
To study the particle number scaling of this superradiant transition and to obtain some analytical insights, we first focus on a mean-field model.The equations of motion for the expectation values s x,y,z ≡ ⟨σ x,y,z ⟩ can be obtained from Eq. ( 3) via the expression ∂ t ⟨ Ô⟩ = Tr([∂ t ρ(t)] Ô), and by performing a mean-field approximation for the two-point correlators ⟨ Â B⟩ = ⟨ Â⟩⟨ B⟩.This results in a set of 3N coupled differential equations governing the dynamics of the form for i and k ∈ {1...N } and pump perpendicular to the x-y plane, such that the atom-dependent driving phase in Eq. ( 1) vanishes (k • r i = 0 for i = 1, ..., N ) (see Appendix D for the more general set of equations).These equations capture the transition between the paramagnetic and magnetized phase remarkably well.In particular, the threshold obtained from the mean-field equations matches with that obtained with the full master equation (3), which are respectively shown by the dash-dotted green trace and the solid black trace in Fig. 2(a).Performing an additional approximation allows us to obtain analytical insights based on this mean-field model.For infinite periodic arrays, we can leverage the geometries' symmetry and define the effective interaction strength and decay rate as , respectively.For subwavelength lattices these effective couplings converge to finite constants as N → ∞, for which the approximation becomes exact.This simplifies the set of equations to just three equations , which describe the mean-field dynamics of a single spin surrounded by all other atoms.They are given as Linearizing this set of equations around the steady-state solution via s x,y,z = s x,y,z 0 + δs x,y,z (t) and solving for the steady-state in δs x,y,z (t) (see Appendix C) allows the extraction of a concise expression for the critical driving strength.In accordance with the Dicke limit [see Fig. 3(a)], we again define the value of Ω where s y is maximal as the critical driving strength Ω c .We obtain (see Appendix C) Despite the substantial approximation performed with the linearization of equations ( 5), this analytical treatment captures the overall properties of the phase diagram, as can be seen from the yellow dashed line in Fig. 2(a).In particular, Ω c has extrema at lattice , the critical driving strength scales linearly with particle number in the Dicke limit, but it doesn't scale with N for the free space case.This feature is independent of lattice spacing, as evinced by the solid grey line (a = 0.2λ0) and the orange dashed line (a = 0.4λ0).
spacings where J eff = 0 (see Fig. 2).These extrema coorespond to a maximum (minimum) when Γ eff > 0 (Γ eff < 0).The critical driving strength given in Eq. ( 6) is constant in the particle number N .This marks a strong difference to the Dicke case, where the critical driving strength in the thermodynamic limit, Ω Dicke c = N γ/2 (see Appendix B), scales linearly with particle number.This distinction between the free space and the Dicke case is illustrated in Fig. 3.The full quantum solution of the Dicke model approaches the analytical steady-state solution (see Appendix B) for N → ∞.For the free space case, we rely on the mean field model in Eq. ( 4) to obtain the particle number scaling.Intriguingly, the critical value Ω c does not scale with particle number, except for finite size effects at very small system sizes[see Fig. 3(b)].The different scaling of the critical driving strength in particle number suggests that the two cases studied here are part of different universality classes.While a detailed analysis of the models' criticality in the thermodynamic limit warrants further study, it goes beyond the scope of the present work.

IV. EMISSION PROPERTIES
One of the core properties of steady-state superradiance in the Dicke limit is the quadratic scaling of the total emission rate γ tot = ij Γ i,j ⟨σ + i σ− j ⟩ in the steadystate (γ ss tot ) with particle number N [34,37].This rate is given as the expectation value of either the jump term or the anti-commuting part of the Lindbladian in Eq. ( 2) via γ tot = Tr where we defined the dissipative Hamiltonian Ĥdis .In Fig. 4(a) we plot γ ss tot obtained from a solution of the master equation in the Dicke limit as a function of N for different driving strengths Ω.For sufficiently large Ω, the total emission rate scales quadratically with particle number.For small Ω, the emission rate saturates to a constant value above a certain particle number [see Fig. 4(a)].This occurs because the decay rates of the states in the symmetric Dicke ladder scale at least with N, and small values of Ω are not enough to sufficiently invert large systems to attain the N 2 scaling characteristic of the |S = N/2, m = 0⟩ state.For subwavelength arrays in free space, where the underlying geometry results in cooperative shifts and decay rates, a natural question arises: How does the geometry influence the emission properties, in particular in the superradiant regime?
Since solving the full master equation ( 3) is not feasible in the free space case, we perform a cumulant expansion up to second order [45][46][47][48] to extract the scaling of the total steady-state emission rate.The full set of equations Similarly, the driving field only brings population back into the fully excited state |E⟩ via the symmetric Dicke ladder, i. e., via state |B⟩.As a result, the steady-state at large drive contains an equal population in |G⟩, |B⟩ and |E⟩.In the free space case (bottom), decay from |E⟩ to |G⟩ can occur either through the bright state |B⟩ at a rate ΓB = γ0 + Γ12, or through the dark state |D⟩ at a rate ΓD = γ0 − Γ12 Consequently, the steady-state exhibits equal population in all four states |G⟩, |B⟩, |E⟩ and |D⟩, which results in a diminished total emission rate.This occurs even in the limit a/λ0 → 0, when Γ12 → γ0 and ΓD → 0. Notably, the intuition gained from the two atom case also holds for larger particle numbers.(b) Overlap of the steady-state with the full Lindbladian spectrum for a chain with lattice spacing a = 0.2λ0 and strong pump Ω = 40γ0 for multiple particle numbers N .In all cases, all states of the spectrum are equally populated in the steady-state.(c) Decay rate κn = −Re(λn) of the eigenmodes of the Liouvillian in Eq. (3) as a function of spacing a.The eigenvector associated with the eigenvalue κn = 0 corresponds to the steady-state.The timescale needed to reach this steady-state is equal to the inverse of the smallest non-zero decay rate (orange markers) of the Liouvillian spectrum.For a/λ0 → 0, this timescale diverges.(d) Emission rate as a function of time for two atoms separate by different distances a, illuminated by a plane wave drive perpendicular to the atomic chain.As a decreases, the time required to reach the steady-state increases and the evolution resembles that of the Dicke limit for longer times.(e) Total emission rate measured at different finite times as a function of the angle of the impinging driving field with the atomic chain.Larger emission rates are attained for perpendicular drive, θ = π/2, than for drive parallel to the atom chain, θ = 0.A spacing of a = 0.1λ0 is considered.can be found in Appendix D. In contrast to the quadratic scaling in the Dicke limit, we find linear scaling with particle number for any value of Ω in the free space case [see Fig. 4(b)], independent of spacing or geometry.This drastic change in the emission properties arises from the space dependent coherent J ij and dissipative Γ ij dipole interactions, which result in additional decay channels that couple bright to subradiant states through dissipation.These decay channels and the subsequent coupling to subradiant states are suppressed in the Dicke limit, which is restricted to the symmetric bright or radiating states [see Fig. 5(a)].While the data presented in Fig. 4 are for a chain of atoms, qualitatively similar results are obtained for a square array or other higher dimensional geometries.
The general mechanism resulting in this stark deviation from the Dicke limit can be already illustrated in a two atom model.For two atoms at a distance a < λ 0 , the dipole-dipole interactions become significant and give rise to coherent interactions J the antisymmetric dark state couples to these states via the suppressed decay rate Γ D = γ 0 − Γ 12 .In the Dicke limit, the distance between emitters tends to zero and the dissipative interaction approaches the spontaneous decay rate, lim |r12|→0 (Γ 12 ) = γ 0 [22,49].As a result, the dynamics is restricted to the symmetric subspace, i. e., the Dicke ladder |E⟩ → |B⟩ → |G⟩ [see Fig. 5(a)], and the population of the dark state is zero at all times (p Dicke D = 0).This changes in the free space case, where an additional dissipation channel to the dark state emerges [50] and Γ D takes a small but finite value, Γ D ̸ = γ 0 .Then, the decay into the dark state is no longer fully suppressed [see Fig. 5(a)], and a finite population in the dark state (p f.s.D ̸ = 0) is attained.For this simple two atom model, the populations of the individual states in the steady-state can be obtained analytically (see Appendix E).In the Dicke case and for sufficiently strong driving on resonance with the bright state, i. e., ∆ = J 12 , the interplay of continuous drive and collective dissipation along the Dicke ladder results in an equillibrium configuration where all states in the symmetric subspace are equally populated, i. e., p Dicke In the free space case, however, the additional dissipation channel ∝ Γ D modifies the equilibrium state, such that all four eigenstates are populated equally in the steadystate, i. e., p f.s.G = p f.s B = p f.s D = p f.s E = 1/4.That is, a significant amount of the population is then trapped in a non-radiative dark state.This results in a reduced emission rate in the free space case, and will ultimately lead to the different scaling of the emission properties with particle number shown in Fig. 4.
Crucially, the steady-state populations at large drive follow the same trend for a general particle number N: in the free space case, all states are equally populated; in the Dicke case, only the states within the symmetric subspace are equally populated.While a full analytical solution for the steady-state of Eq. ( 3) is cumbersome in this general setting, we can numerically test this intuition for small atom numbers.In particular, the steady-state solution can be determined via the spectrum , where κ n and ν n respectively denote the decay rate and energy shift associated to the n-th eigenvalue.More precisely, the steady-state fulfills ∂ t ρss = 0 and therefore corresponds to the Liouvillian eigenstate with zero decay rate, i. e., the state with κ n = 0.In Fig. 5(b), we show the overlap O L = ⟨ψ| ρss |ψ⟩ of the steady-states with the individual eigenstates |ψ⟩ of the many-body Hamiltonian Ĥ for the free space case.As occurs in the two particle case, the steady-state has equal overlap with all 2 N eigenstates of the Hamiltonian for all simulatable particle numbers up to N = 8.The steady-state emission rate, can be expressed as the expectation value of the dissipative Hamiltonian Ĥdis = i,j Γ ij σ+ i σ− j .Then, each of the 2 N eigenstates of H dis can be assigned to one of the N +1 excitation subspaces containing m ∈ {0, N } excitations.There are N m = M such states in the m-excitation subspace, each having a decay rate Γ (m) i , and such that the sum of all decay rates within the m-th excitation subspace is equal to Since all eigenstates in all excitation manifolds are equally populated, the total emission rate is simply given as the average of all their decay rates (7) This confirms the linear scaling of the total emission rate with particle number obtained in Fig. 4(b).Intuitively, the fact that all eigenstates of H dis are equally populated results in significant contribution of dark decay rates with Γ D ≪ γ 0 in the sum of Eq. ( 7), which drastically diminishes the total emission rate.
In contrast, only the N + 1 states contained in the symmetric subspace are occupied in the Dicke limit.The decay rate of the symmetric state in the m-excitation subspace, |S = N/2, m⟩, is given as , and the total emission rate then reads which highlights the quadratic scaling in particle number in the Dicke case.These findings are in agreement with the numerical results shown in Figs.4(a) and (b).Hence, the stark difference in emission properties between free space and the paradigmatic Dicke limit can be explained by the contribution of the different dissipation channels available in the two cases.While only radiative or superradiant states contribute to emission in the Dicke case, non-radiative dark or subradiant states are significantly occupied in the free space case, resulting in diminished photon emission.

V. FINITE TIME EFFECTS
The emergence of decay channels with non-zero decay rates much smaller than γ 0 affect not only the steadystate emission properties of the system, as discussed in Sec.IV, but also its dynamics.So far, we have characterized the steady-state properties by either determining the null-space of the Liouvillian or evolving approximate equations of motion for very long times (t > 100γ 0 ).In experiments, however, the measurement of the emission properties typically takes place at much earlier times, and it is important to understand the interplay between the measurement time and the slowest or characteristic timescale at which the system evolves.To do so, we note that the time evolution of the density matrix ρ under the master equation ( 3) can be expressed as ρ(t) = 2 N n=1 c n e λnt u n , where the coefficients c n are fixed via the initial condition, and λ n and u n respectively denote the eigenvalues and eigenvectors of the Liouvillian.Note that the eigenvalues are typically complex, λ n = −κ n + iν n , and are characterized by a decay rate κ n and an energy shift ν n .Then, the fundamental timescale at which the steady-state is reached, τ ss , is equal to the inverse of the smallest non-zero decay rate of the Liouvillian spectrum with a non-zero contribution c n .In Fig. 5(c), we show the decay rates of the Lindbladian spectrum for two atoms, κ n = −Re(λ n ).The smallest non-zero decay rate is indicated by yellow markers.For decreasing lattice spacing a → 0, this eigenvalue approaches zero, which implies a divergence of the timescale required to reach the steady-state.This effect is nicely illustrated in Fig. 5(c), where we show the time evolution of the total emission rate for two atoms for different atom distances.For large enough lattice spacings, the steady-state emission rate reaches the analytic free space value given in Eq. ( 7) very quickly.For decreasing lattice spacings, the emission rate at finite times gets closer and closer to the Dicke value given in Eq. ( 8), as it takes longer and longer times for the subradiant states to be populated.For a lattice spacing of a = 0.05λ 0 and after some oscillatory initial dynamics, the time evolution of the emission rate overlaps with that of the Dicke limit for the time window shown, t ∈ {0, 20γ 0 }.For such small lattice spacings, the subradiant decay rate is heavily suppressed and the time τ ss required to reach the actual free space steady-state emission rate becomes increasingly large.
In other words, the emergence of decay channels with heavily suppressed but non-zero decay rates increases the time required for the system to equilibrate and populate its dark states.This modifies the emission properties when measured at finite times (as it would occur in any experimental implementation).In Fig. 6, we illustrate the effect of finite-time measurements on the total emission rate.As shown in Fig. 6(a), the emission rate approaches that of the Dicke limit for smaller interparticle spacings (i.e., larger densities) and earlier measurement times.
Another important parameter which we did not consider thus far is the angle of the incident pump beam with respect to the orientation of the lattice.While the final steady-state is independent of the angle of incidence, the transient dynamics leading to this steady-state strongly depends on the phase factors ∝ e ik•ri in Eq. (1).For the two atom case, we solve the full master equation [see Fig. 5(e)] and find that the effect of finite measurement times is most pronounced for perpendicular illumination (θ = π/2), where the total emission rate gets close to the actual Dicke value at early times.For parallel illumination (θ = 0.0), this effect is less pronounced, and we obtain emission rates much closer to the free space case for all measurement times.This phenomenon can be understood by noting the fundamental difference between both scenarios.For perpendicular illumination, the driving field couples the ground state |G⟩ only to  For a measurement at finite times, the values for γtot approach the Dicke limit for decreasing density.(b) Scaling of the emission rate with particle number for a fixed chain length L = 3λ0 and varying particle number, for different measurement times t0 and plane wave drive perpendicular to the chain axis.The linear black dashed line correspond to measurement times γ0t = 100, for which the emission rate is consistent with the the steady-state analysis shown in Fig. 4. For early measurement times, the scaling becomes superlinear and therefore resembles the Dicke case show in Fig. 4(a).(c) Same as (b), but for plane-wave drive parallel to the chain axis.Compared to perpendicular drive, the emitted intensity is closer to the free space case and its characteristic linear scaling for all measurement times.(b) and (c) therefore show that the properties obtained for the two atom model in Fig. 5(e) also hold for larger particle numbers.
the bright state |B⟩, and contributes to a rapid population of the superradiant states of the system.The subradiant state |D⟩, however, only gets populated via the subradiant decay channel at a rate proportional to the slow timescale Γ D = γ 0 − Γ 12 .For parallel illumination, however, the driving field generally couples the ground state |G⟩ to both the superradiant |B⟩ and subradiant |D⟩ single-excitation states, and the population of the latter can occur at a much faster timescale.As a result, the steady-state is also reached at earlier times, and the effect of measuring at finite times is substantially suppressed.
We further study the scaling of the emission rate as a function of density.To do so, we fix the length of the chain and analyze the total emission rate as a function of particle number N for different measurement times [see Fig. 6(b)-(c)].We again simulate the dynamics of the system by means of a second order cumulant expansion (see Appendix D) given the large particle numbers considered.For sufficiently late measurement times (γ 0 t 0 = 100), we retrieve the expected linear scaling char-acteristics of the steady-state analyzed in Sec.IV.For earlier measurement times, however, the total emission rate scales superlinearly with particle number, confirming that the effect observed in the minimal two atom model also occurs for large system sizes.Again, the direction of the plane-wave driving field has the same effect as in the two-atom case.The strongest deviation from linear scaling is found for perpendicular drive, which predominantly excites the superradiant modes of the system.For driving along the chain axis, the coupling to the most superradiant modes is reduced [51] and deviations from linear scaling can only be observed at very short measurement times.

VI. CONCLUSIONS
Motivated by the recent experimental progress in generating controllable subwavelength geometries of quantum emitters, we analyzed the role of geometry induced shifts and decay rates on the steady-state properties of driven atom arrays.We particularly focused on the transition between a magnetic and a paramagnetic (superradiant) regime, and characterized the emission properties for varying geometries.We find that the emergence of collective shifts and dissipations in structured arrays results in a stark difference compared to the simple Dicke regime, where equal all-to-all interactions with vanishing shifts are assumed.In particular, the population of dark states destroys the quadratic scaling with particle number of the steady-state emitted light intensity, characteristic of the Dicke limit.Instead, only a linear scaling is observed for extended geometries.At finite measurement times and for very dense arrangements of atoms, however, we approach the Dicke-like behavior due to the diverging time scales to reach the actual steady-state.These insights could directly contribute to the interpretation of recent experimental observations in a dense driven pencilshaped cloud of atoms, which experimentally realized the driven superradiant phase transition in free space [34].
In a broader setting, driven dissipative spin models are an exciting avenue for future research.The control over various atomic degrees of freedom enables novel protocols and phenomena based on the interplay between dissipation and coherent drive.Such protocols could allow the realization of novel states of matter in the steadystate, and provide alternative routes for dissipation assisted state preparation of light and matter.The latter could be achievable in generalized analog settings as considered in this work, as well as in fully programmable quantum devices, where Lindbladian terms can in principle be engineered by performing appropriate randomized measurements to target a particular steady-state.
Calculating additional equations of motion for the two-point correlators ∝ ⟨O i O k ⟩ and replacing averages over third-order operators by [45,46]

Figure 2 .
Figure 2. (a) Steady-state phase diagram for a square with the collective spin pointing in z-direction, ⟨Sz⟩/N , as a function of Rabi frequency Ω and lattice spacing a for a four atom square.The solid black line indicates the numerically obtained threshold from the full quantum solution, the green dash dotted line indicates the threshold obtained via the mean-field equations, and the dashed yellow line shows the analytically obtained threshold given by Eq. (6).The horizontal dashed white lines indicate the lattice spacings at which the effective coupling strengths J eff shown in (b) is zero.(b) Effective coupling strength J eff and decay rate Γ eff as a function of lattice spacing.The vertical dashed lines indicate zero crossings of J eff .

Figure 3 .
Figure 3. (a) steady-state value of the collective spin expectation value ⟨ Ŝα⟩ (α ∈ {x, y, z}) obtained by solving the master equation (3) for the steady-state in the Dicke limit as a function of driving strength Ω.The thick black lines show the analytical solution obtained from a mean-field model (see Appendix B).(b) steady-state solutions of the mean-field model given in Eq. (4) for a square lattice in free space with lattice spacing a = 0.2λ0.The different shadings in (a) and (b) indicate different particle numbers N = 4, 36, 64, 100 from light to dark.In either case, we use the maximum of ⟨Sy⟩ as the threshold condition.(c) Scaling of the critical driving strength Ωc as a function of particle number.As suggested by (a) and (b), the critical driving strength scales linearly with particle number in the Dicke limit, but it doesn't scale with N for the free space case.This feature is independent of lattice spacing, as evinced by the solid grey line (a = 0.2λ0) and the orange dashed line (a = 0.4λ0).

Figure 4 .
Figure 4. Scaling of the total emission rate as a function of particle number (a) for the Dicke case, and (b) for a chain with lattice spacing a = 0.2λ0 in free space, obtained by evolving the equations in second order cumulant expansion until a steady-state is reached.The different colors and markers indicate different driving strengths.For sufficiently strong driving, we find a quadratic scaling ∝ N 2 with particle number for the Dicke case and a linear scaling ∝ N for the free space case.

Figure 5 .
Figure 5. (a) Intuitive explanation for the different scaling of the total emission rate with the particle number N based on the two atom case.In the Dicke case (top), decay from |E⟩ to |G⟩ only occurs via the bright state |B⟩ at a rate ΓB = 2γ0.Similarly, the driving field only brings population back into the fully excited state |E⟩ via the symmetric Dicke ladder, i. e., via state |B⟩.As a result, the steady-state at large drive contains an equal population in |G⟩, |B⟩ and |E⟩.In the free space case (bottom), decay from |E⟩ to |G⟩ can occur either through the bright state |B⟩ at a rate ΓB = γ0 + Γ12, or through the dark state |D⟩ at a rate ΓD = γ0 − Γ12 Consequently, the steady-state exhibits equal population in all four states |G⟩, |B⟩, |E⟩ and |D⟩, which results in a diminished total emission rate.This occurs even in the limit a/λ0 → 0, when Γ12 → γ0 and ΓD → 0. Notably, the intuition gained from the two atom case also holds for larger particle numbers.(b) Overlap of the steady-state with the full Lindbladian spectrum for a chain with lattice spacing a = 0.2λ0 and strong pump Ω = 40γ0 for multiple particle numbers N .In all cases, all states of the spectrum are equally populated in the steady-state.(c) Decay rate κn = −Re(λn) of the eigenmodes of the Liouvillian in Eq. (3) as a function of spacing a.The eigenvector associated with the eigenvalue κn = 0 corresponds to the steady-state.The timescale needed to reach this steady-state is equal to the inverse of the smallest non-zero decay rate (orange markers) of the Liouvillian spectrum.For a/λ0 → 0, this timescale diverges.(d) Emission rate as a function of time for two atoms separate by different distances a, illuminated by a plane wave drive perpendicular to the atomic chain.As a decreases, the time required to reach the steady-state increases and the evolution resembles that of the Dicke limit for longer times.(e) Total emission rate measured at different finite times as a function of the angle of the impinging driving field with the atomic chain.Larger emission rates are attained for perpendicular drive, θ = π/2, than for drive parallel to the atom chain, θ = 0.A spacing of a = 0.1λ0 is considered.

Figure 6 .
Figure 6.(a) Total emission rate as a function of inverse lattice spacing for two atoms.For a measurement at finite times, the values for γtot approach the Dicke limit for decreasing density.(b) Scaling of the emission rate with particle number for a fixed chain length L = 3λ0 and varying particle number, for different measurement times t0 and plane wave drive perpendicular to the chain axis.The linear black dashed line correspond to measurement times γ0t = 100, for which the emission rate is consistent with the the steady-state analysis shown in Fig.4.For early measurement times, the scaling becomes superlinear and therefore resembles the Dicke case show in Fig.4(a).(c) Same as (b), but for plane-wave drive parallel to the chain axis.Compared to perpendicular drive, the emitted intensity is closer to the free space case and its characteristic linear scaling for all measurement times.(b) and (c) therefore show that the properties obtained for the two atom model in Fig.5(e) also hold for larger particle numbers.