Signatures of optical phase transitions in super- and subradiant arrays of atoms

Resonant light interacting with matter can support different phases of a polarizable medium, and optical bistability where two such phases coexist. Here we identify signatures of optical phase transitions and optical bistability mapped onto scattered light in planar arrays of cold atoms. Methods on how to explore such systems in superradiant, and extreme subradiant states existing outside the light cone, are proposed. The cooperativity threshold and intensity regimes for the intrinsic optical bistability, supported by resonant dipole-dipole interactions alone, are derived in several cases of interest analytically. Subradiant states require lower intensities, but stronger cooperativity for the existence of non-trivial phases than superradiant states. The transmitted light reveals the onset of phase transitions and bistability that are predicted by mean-field theory as large jumps in coherent and incoherent signals and hysteresis. In the quantum solution, traces of phase transitions are identified in enhanced quantum fluctuations of excited level populations.


I. INTRODUCTION
Resonant emitters in regular planar arrays have attracted considerable attention from classical circuit resonators forming metamaterials [1] and metasurfaces [2] to plasmonics [3] and quantum systems, such as superconducting SQUID rings [4] and cold atoms [5]. Such surfaces can be utilised for manipulation of electromagnetic fields, including phase-holography [6] and sensing [7]. In systems where the radiative interactions between closelyspaced emitters are particularly strong, the entire array has been driven to a giant subradiant state [5,8].
Most experiments on collective optical responses of cold atoms so far, in both random atomic ensembles [24][25][26][27][28][29][30][31][32][33][34][35] and in arrays [5] and chains [36] with unit occupancy, have focused on the limit of low light intensity (LLI), where the full quantum model can, under appropriate conditions, be reduced to a linear system of N harmonic oscillators [37]. Beyond the LLI regime with multiple excitations, atomic arrays start experiencing saturation, and the rich phenomenology of long-range interactions and collective behaviour can lead to the full many-body quantum solutions deviating from the semiclassical models that neglect quantum fluctuations [20]. The differences between quantum and classical solutions in nonlinear systems are widely studied in the context of phase transitions, and in optics, one of the best-known phase transitions is optical bistability [38] in atomic systems.
Optical bistability and phase transitions have been actively studied in systems without the spatial correlations and structure of the sample [39][40][41][42][43][44][45][46][47], e.g., in cavities where the feedback mechanism is provided by the cavity mirrors [48][49][50]. Intrinsic bistability is a process where phase transitions are generated by the self-interactions of the sample, and despite having been observed in highlyexcited Rydberg atoms in the microwave regime [51], intrinsic bistability was for a long time considered unachievable for atoms with light-mediated interactions. Recent theoretical studies that also take into account the spatial structure of the many-body systems suggest that intrinsic bistability and phase transitions are more generic and could occur in a variety of systems with shortand long-range interactions [52][53][54].
For optical systems, it is natural to ask what are the observable signatures of phase transitions and optical bistability, and how these are mapped onto the scattered light. In this paper, by studying light emission from radiatively strongly coupled atoms in planar arrays of subwavelength spacing, we identify optical signatures of phase transitions in collective atomic excitations. To do so, we employ periodic boundary conditions and a meanfield approximation which closes the spectral gap in the system, representing a decohered quantum state where the correlations are absent, and compare our results to the full quantum model.
We develop a simple analytic theory for an intrinsic optical bistability due to radiative interactions between atoms in planar arrays and derive the cooperativity parameter, indicating a bistability threshold ka < (π/3) 1/2 , with the lattice spacing a and resonance wavenumber k. We find that multiple mean-field-theoretical stable phases, including ones with spontaneous symmetry breaking and persistent oscillations, and optical bistability are identifiable in the transmitted light as large jumps in coherent and incoherent signals and hysteresis upon sweeping of the laser frequency. If the corresponding changes in dipole amplitudes are small, the signal of phase transitions and hysteresis in the coherent transmission is sometimes much weaker than in the incoherent photon count, which still provides sharp peaks, e.g., when moving into regions of antiferromagnetic and oscillatory phases. In the quantum solution, the phase transitions and bistability are absent, but traces of them are revealed in enhanced quantum fluctuations of incoherently scattered light, and most clearly in those of the excited level populations.
We find that the response sensitively depends on the underlying LLI collective excitation eigenmode of the corresponding linear system that is targeted by incident light. Bistability can even exist between subradiant and superradiant modes, providing a method for also preparing subradiant excitations via a laser frequency sweep. For subradiant modes, bistabilities occur at lower intensities and the existence of phase transitions requires smaller lattice spacings, ka 0.34π, compared to the one for a superradiant mode, ka 0.44π. We propose methods on how to drive such eigenmodes by manipulating the atomic level shifts and consider two examples: a uniform mode that was recently experimentally studied in subradiant transmission measurements [5], which at smaller lattice spacings, considered here, becomes superradiant, and an extreme subradiant checkerboard eigenmode that can exist outside the light cone, decoupled from the environment.

A. Quantum system of atoms and light
We consider a two-level system of cold atoms trapped in a two-dimensional (2D) array with one atom per site, illuminated by an incident plane wave E + (r) = E 0ê exp(ikz); Fig. 1. We take the polarisation and the direction of the atomic dipoles to beê = −(x +ŷ)/ √ 2 along the diagonal of the lattice. Light-induced resonant dipole-dipole interactions mediate strong interactions between the atoms. The atomic systems are subject to periodic boundary conditions to simulate an infinite lattice and, for simplicity, we vary the parameters of only N = 4 atoms in a square array. We also assume that the atoms are sufficiently tightly confined, such that the spatial fluctuations can be neglected. The standard many-body quantum master equation for the atoms in the rotating-wave approximation for slowly varying amplitudes reads [55] where the square brackets represent a commutator and σ + j = |e j j g| = (σ − j ) † the raising operator, where |e j and |g j are the excited and ground state of the two-level atom on site j, respectively. The dispersive and dissipative parts of the light-induced dipole-dipole interaction terms are Ω jl and γ jl , respectively (see Appendix A for

details). The Hamiltonian is given bŷ
eg is the detuning, ω = kc the laser frequency, ω (l) eg the transition frequency of an atom on site l, and d eg the dipole matrix element, with d ge = d * eg . We express the incident light intensity I = 2 0 c|E 0 | 2 in units of the saturation intensity, I sat = c4π 2 γ/3λ 3 , or the Rabi frequency, R l = d eg · E + (r l )/ , as I/I sat = 2(R/γ) 2 , where the single-atom linewidth γ = |d eg | 2 k 3 /(6π 0 ).

B. Mean-field approximation
In addition to the full quantum many-body dynamics, we also consider the Gutzwiller mean-field approximation,ρ ≈ ⊗ρ i , where quantum fluctuations between the atoms are neglected. This corresponds to the factorization of internal level correlations, σ α iσ β j ≈ σ α i σ β j (α = β), since we assume atoms are at fixed positions with no spatial fluctuations, and therefore there are no light-induced correlations [20,37] between the atoms after the factorization. The dynamics then obey the non-linear equationṡ where ρ (l) ge = Tr{σ − lρ (t)} and ρ (l) ee = Tr{σ ee lρ (t)}. We will use Eqs. (3) to determine the long-time phases and optical bistability that can occur in the system.

C. Scattered light
The total light amplitude is the sum of the incident and scattered fieldsÊ ± (r) = E ± (r) +Ê ± s (r), with the scattered electric field given by the sum of the contributions from all the atoms 0Ê + where G(r−r l ) is the dipole radiation kernel [56] [Eq. (13) in Appendix A]. We will compare the optical responses obtained from the full quantum dynamics of Eq. (1) with those calculated from the mean-field equations, Eqs. (3). We consider coherently transmitted light in the forward direction, T coh = |ê · Ê − (r) | 2 / |ê · E − (r)| 2 , expressed in terms of the optical depth OD = − log(T coh ). We also calculate the rate of scattered photons where Ê − s (r) ·Ê + s (r) = Ê − s (r) · Ê + s (r) for the coherent and Ê − s (r)·Ê + s (r) − Ê − s (r) · Ê + s (r) for the incoherent photon count-rate (ICR) [Eq. (15) in Appendix A]. The ICR that we will use under the mean-field description [Eq. (16) in Appendix A] is different from the usual semiclassical approximation for the incoherent scattering. The atom-light dynamics is solved from the meanfield equations, Eqs. (3), but the single-atom quantum description of emitted light σ ee l is now included [20] for the scattered light, and the ICR no longer vanishes for atoms at fixed positions [57].

III. OPTICAL SIGNATURES OF PHASE TRANSITIONS AND BISTABILITY
A. Collective low light intensity eigenmodes In the limit of LLI, ρ (l) ee = 0, and the mean-field equations, Eqs. (3), coincide with the coupled-dipole model of classical linear oscillators driven by light. In this regime we may analyse the optical response using LLI collective radiative excitation eigenmodes and the complex eigenvalues, which represent the collective line shifts (from the resonance of an isolated atom) δ q and linewidths υ q . Collective modes with υ q > γ (υ q < γ) are termed superradiant (subradiant). Modes with |q|a > 0.2π exist outside the light cone and are completely dark in an infinite lattice.
We focus on two LLI eigenmodes (see Appendix B): the spatially uniform superradiant mode v un (r l ) (υ uni = 25γ) and a subradiant mode with a checkerboard-patterned phase variation with every atom oscillating π out-ofphase from its nearest-neighbour v cb (r l ) (υ cb = 2 × 10 −4 γ). To simulate an infinite system, we use periodic boundary conditions by adding repeat images of the system to the boundaries. We truncate to 101 images along thex andŷ direction, which gives an effective lattice size of 406×406, and non-zero linewidth for the checkerboard mode.
The uniform eigenmode v un directly couples to the normally incident light of uniform phase profile, resulting in a broad resonance in the OD, shown in Fig. 2(a). We show that we can also excite the checkerboard eigenmode v cb and prepare coherent strongly subradiant excitations that exist outside the light cone. This can be achieved by ac Stark shifts [58] of lasers (or microwaves) forming a checkerboard pattern of atomic level shifts that results from a standing-wave, cos 2 [5k(x+y)/ √ 2], with the intensity varying along the lattice diagonalx +ŷ and the intensity maxima separated by √ 2a. Alternating blue-and red-detuned atomic transitions for adjacent atoms cause them to oscillate π out-of-phase, resulting in the excitation of the checkerboard subradiant eigenmode. The relative angle between the field generating the ac Stark shift and the lattice can be adjusted to control the periodicity. An example of an atomic transition particularly suitable for closely-spaced atoms is 3 P 0 → 3 D 1 in 88 Sr [59], which can have a resonance wavelength of λ 2.6µm and spacing of 206.4nm, resulting in the effective lattice spacing a 0.08λ. Figure 2(a) shows how these checkerboard-patterned alternating level shifts lead to a coupling to the checkerboard subradiant mode, producing a Fano resonance in the OD at ∆ = 10.8γ, with the corresponding large population of the checkerboard eigenmode at this resonance [ Fig. 2

B. Analytic results for optical bistability
Classifying the steady states of the mean-field solutions of Eqs. (3) determines the phases that emerge as a function of detuning and incident intensity. We calculate the general phase diagram numerically. However, it is important to understand the collective effects in optical bistability by first deriving solutions in some special cases analytically. In order to do so, we consider the uniform case by substituting ρ  into Eqs. (3). We then obtain the stationary states where we have defined which, withΩ = j =l Ω jl andγ = j =l γ jl , is the total external electric field (incident plus scattered field from all the other atoms, given in terms of the Rabi frequency) driving an arbitrary atom l in the ensemble. Eqs. (6) are equivalent to the familiar solutions of the independentatom optical Bloch equations [Eqs. (28) in Appendix C], but with the Rabi frequency R replaced by R eff . As ρ ge appears on the both sides of Eq. (6a) via R eff , we generally have multiple solutions. For two different coexisting stable solutions, we have optical bistability. We can eliminate from Eqs. (6a) and Eq. (7) the atomic variables and obtain an equation for the incident light field R = R(R eff ). The bistability threshold is then found when d|R| 2 /d|R eff | 2 = 0. This gives a cubic polynomial for |R eff | 2 [see Eq. (30) in Appendix D] in terms of γ, ∆,Ω, andγ. Simple analytic expressions for the optical bistability threshold can then be obtained for ∆/γ =Ω/γ, yieldingγ > 8γ, and for ∆/γ = −γ/Ω, yieldingΩ 2 > 27γ 2 . Below these values, there is no bistability for any intensity. We can also obtain analytic forms for the bistable solutions of the external field R eff acting on an atom, forΩ,γ ∆ 2 , where we have defined the single-atom excited state occupation for unsaturated drive, p unsat = |R| 2 /(∆ 2 + γ 2 ), and the cooperativity parameter, The two solutions represent very different responses to the incident light. The first "cooperative" solution, Eq. (8a) (in an analogy with the terminology of optical bistability in cavities [41]), exists for p unsat < |C| 2 /2 and arises due to the atoms behaving collectively, creating a field that counteracts the incident light and resulting in the atoms absorbing strongly, with enhanced absorption for larger atom density. The second "single-particle" solution, Eq. (8b), exists when p unsat > 2(|C| + Re[C]), and arises when the atoms react to the incident light almost independently, with R ≈ R eff when |R| → ∞.
The atoms now saturate and absorption is weak, with the medium becoming transparent.
The simplest system exhibiting collective interactions is that of two atoms (Ω = Ω 12 ,γ = γ 12 ). In this case, we can satisfyΩ 2 > 27γ 2 for closely spaced atoms for ∆/γ = −γ/Ω. Approximating the resonant dipole-dipole coupling by Ω 12 ∼ 1/(ka) 3 , where a denotes the atom separation, results in the bistability threshold of roughly ka 1, with the precise value depending on the orientation of the dipoles (see Appendix E).
Analytic expressions can be obtained for atomic chains and arrays for ∆/γ =Ω/γ, where the bistability threshold is independent ofΩ. For an infinite 1D chain, we can sum the series of dissipative dipole-dipole interaction terms over the atoms to obtain the collective resonance linewidthγ wherer indicates the atomic chain orientation. The bistability thresholdγ > 8γ is met when ka < π/6 (a 0.08λ) or ka < π/12 (a 0.04λ) for dipoles parallel and perpendicular to the chain, respectively. For an infinite 2D array, with a uniform distribution of atomic dipoles in the plane, we obtain [60] (see also Ref. [61]) for the collective linewidthγ 2D /γ = 3π/(ka) 2 − 1, which allows for larger lattice spacings, ka < (π/3) 1/2 (a 0.16λ), for the bistability threshold than a 1D chain. For dipoles normal to the plane,γ = −γ, and soγ > 8γ and bistability is not possible. The analogy between the optical bistability in atom arrays and in cavities [39][40][41][42][43][45][46][47] can now be most easily illustrated, and our adapted terminology motivated, at the specific value of ∆/γ =Ω/γ for which C =γ/2γ in Eq. (9) is real. The expression for the incident light field [Eq. (29) in Appendix D] then has a similar form as that in cavity systems, with the same formulaic dependence on the cooperativity parameter in atom arrays as that for optical bistability in cavities [41]. In cavity QED, the cooperativity parameter represents recurrent interactions of an atom with light reflecting between the cavity mirrors. In atom arrays for ∆/γ =Ω/γ, the bistability condition C 4 then translates to the density threshold ka ∼ 1 -equivalent to the requirement for the existence of substantial recurrent and correlated light scattering, where the light is scattered more than once by the same atom [37,62,63]. Moreover, asγ in atom arrays takes the role of the atom-cavity coupling coefficient, the condition C 1 then also corresponds to the strong coupling regime of cavity QED.
Numerical solutions of the phase diagram in a planar array agree well with the analytic result of the optical bistability ka < (π/3) 1/2 for the uniform phases when driving the superradiant eigenmode, and for ka 0.44π, only one phase persists and no phase transitions occur. When specifically targeting the subradiant eigenmode by using alternating checkerboard-patterned level shifts, optical bistability can only be predicted by numerically solving the equations of motion [specifically Eqs. (24) in Appendix C], with a much smaller spacing ka 0.34π needed for the optical bistability and phase transitions to occur. Bistability also occurs at lower intensities, with bistability in the range 0.07 I/I sat 190 when driving the subradiant mode compared to 2 I/I sat 406 when driving the superradiant eigenmode.

C. Uniform level shifts: Mean-field phases and optical signatures
So far we have studied the coupling of light in the limit of LLI and the emergence of optical bistability for a uniform atom array. Next we determine the entire phase diagram of atoms coupled by light-mediated interactions beyond the LLI regime by finding the steady-state solutions of Eqs. (3) for ka = 0.2π. In general, we find the system can exhibit spatially uniform phases, antiferromagnetic (AFM) phases and persistent oscillations (OSC), as well as different optical bistabilities. Transitions between different phases can result in small dips in the OD, and lead to large peaks in the ICR. Phase bistabilities can result in large jumps in the OD and ICR, as well as hysteresis upon varying the laser frequency.
We first consider the case where the atomic level shifts are all equal; Fig. 1(a). Beyond the obvious phases representing uniform low and high excitation numbers, labelled U 1 and U 2 [the q = 0 case of Eq. (25) in Appendix C], respectively, we interestingly also find stable phases with spontaneously broken translational symmetries and regions of two coexisting stable phases. The coherent and incoherent optical responses, OD and ICR, from the mean-field analysis are shown in Fig. 4.
While the uniform phases U 1 and U 2 vary smoothly into one another [white regions of Fig. 1(a)], there also exists a U 1 /U 2 bistability due to two possible values of ρ ee , where the state of the system depends on the initial condition (dark blue region). This bistability region is largely well described by our earlier analytics and contours in Fig. 3 [derived from Eq. (29) in Appendix D]. However, differences occur due to one of the uniform phases becoming unstable at positive detunings. The broad resonances of U 1,2 for I/I sat = 100 in Fig. 4 correspond to the underlying superradiant LLI uniform excitation eigenmode [Eq. (20) in Appendix B]. For the U 1 /U 2 bistability at I/I sat = 200, we find hysteretic behaviour upon sweeping the resonance from either red-or blue-detuned side. This demonstrates how crossing a region of bistability results in a large jump in both the OD and ICR, with the jump point depending on the initial condition.
Despite the uniformly excited atoms, stable phases with spontaneously broken translational symmetries emerge with the atomic dipoles oscillating π out of phase in the neighbouring sites [red regions in Fig. 1(a)]. These AFM phases appear at detunings resonant with the LLI excitation eigenmodes u ±,cb [Eqs. (21) and Eq. (22) in Appendix B], and have the same underlying spatial variation, with a striped AFM ± phase originating at ∆ = 4.65γ and a checkerboard AFM cb phase originating at ∆ = 10.8γ. The AFM phases materialise as nonlinear interactions between the atoms allow small fluctuations to populate the spatially nonuniform modes in the system, causing phase instabilities. Two narrow peaks for I/I sat = 100 in Fig. 4 in the ICR at ∆ = 3.8γ and 8.7γ signal spontaneous symmetry breaking and a phase transitioning from U 1 to the AFM ± and AFM cb , respectively. No clear signature of this transition can be seen in the OD. We find that the AFM phase is bistable with the U 1 phase [yellow regions of Fig. 1(a)]. Switching off the incident drive, the uniform phase decays superradiantly, and the AFM phase decays subradiantly. Therefore, AFM cb /U 1 bistability represents an interesting situation where either a superradiant or subradiant phase can be populated depending on the initial condition. The hysteresis associated with the bistability could be utilised as a novel method for preparing subradiant excitations, when a steady-state superradiant mode is transformed into a subradiant one by a laser frequency sweep.
The AFM phases also can become unstable (via Hopf bifurcations) resulting in a phase where the atoms continue to oscillate indefinitely, denoted OSC. Such phases appear as additional peaks at I/I sat = 100 and I/I sat = 200. Due to oscillations in the OSC phase, the signal is noisy as the stationary state is no longer well defined. There are also regions of OSC and U 1 bistability. For the OSC/U 1 bistability at I/I sat = 100, clear hysteresis can be seen in the ICR, but hysteresis in the OD is very small. This is due to the alternating out-of-phase dipoles, which, when summed together to calculate the OD, nearly cancel. The ICR always shows clear peaks and hysteresis as it depends on the excitation strength and not the phase. There are small regions (not marked) near this OSC/U 1 bistable region where a new phase emerges which is neither spatially uniform or AFM in nature, where two dipoles are out of phase to one another, and the two remaining ones in phase with one other.

D. Uniform level shifts: Quantum fluctuations
In the full quantum theory, there is no bistable behaviour or phase transitions. While generally at high intensities the mean-field and quantum results are in closer agreement [20] as the atoms start to scatter more, at intermediate intensities we find in Fig. 4 considerable deviations where the mean-field solutions display bistability. In the full quantum description, due to quantum correlations between different atoms, the ICR no longer represents the excited level population as in the single-atom quantum description that we use to analyse the meanfield dynamics. However, the ICR still shows a resonance near the detunings where mean-field AFM transitions occur.
It is generally known from past bistability studies [44-46, 50, 54, 64] that mean-field bistabilities coincide with enhanced quantum fluctuations, which can be under- are strongly enhanced around the U 1 /U 2 and AFM cb /U 1 phase bistabilities. This corresponds to large variations in the excitation strength between the different meanfield solutions, and could be detected by resonantly transferring excited atoms to another level.

E. Alternating level shifts: Mean-field phases and optical signatures
By engineering a checkerboard pattern of alternating atomic level shifts, detailed earlier in the discussion of the LLI modes, we are able to drive collective excitations where the atoms oscillate π out-of-phase with respect to their nearest-neighbour, and whose LLI limit represents subradiant checkerboard eigenmode [Eq. (22) in Appendix B] existing outside the light cone. We now analyse how this influences the phase diagram beyond the LLI limit. The alternating level shifts explicitly break the translational symmetry of the lattice. The spatially uniform phases U 1,2 of Fig. 1(a) now transform to checkerboard AFM cb phases in Fig. 1(b), but can still be distinguished as having low and high excitations, labelled AFM 1 and AFM 2 . AFM phases no longer occur in pairs as the level shifts favour one AFM configuration over the other. The coherent and incoherent optical responses, OD and ICR, from the mean-field analysis are shown in Fig. 6. One key difference from the uniform level shifts is that the bistability now occurs at lower intensities, as discussed earlier when analysing the bistability analytics. No bistability is found for I/I sat = 200, so instead we look at I/I sat = 180. Deviations from the LLI model and the emergence of nonlinear response depend on the linewidth of the corresponding LLI eigenmode, with subradiant modes being more sensitive at lower intensities than superradiant ones [65]. Therefore, for nonuniform level shifts, the intensities at which non-trivial phases emerge are lower due to the checkerboard subradiant mode [Eq. (22) in Appendix B] being populated.
There is still a region of OSC phase and also a small region of AFM 1 /OSC bistability, indicated by several corresponding peaks in the ICR (and also a dip in the OD for I/I sat = 180). There are a few cases where the OSC  Fig. 1(b) phase becomes unstable and only the AFM 1 phase persists, which are not marked. Subradiant excitations in the limit of LLI lead to narrow Fano resonances when interfering with broader-resonance modes, as shown in Fig. 2. Some of these transform to bistable regions, such as AFM 1 /AFM 2 , which again displays large jumps and hysteresis upon sweeping the detuning.
Finally, there are regions (green) where phases emerge that are not spatially uniform or AFM in nature (labelled BD 1,2 ), and are bistable with an AFM phase. For the BD 1 (BD 2 ) phase, the atoms alongx +ŷ are in-phase (out-of-phase), while the atoms alongx −ŷ are out-ofphase (in-phase). Small regions of the BD 2 phase were found for the uniform shift case. Within the AFM/BD 2 bistability region, the BD 2 phase can become unstable and only the AFM phase remains. BD 1,2 phases occur because of the striped subradiant modes [Eqs. (21) in Appendix B], which have a spatial variation that does not match with the level shift profile, but are populated by nonlinear interactions even though they do not couple to the drive. There is a small peak in the ICR at ∆ = 0.9γ in Fig. 1(b) for the AFM 1 to BD 1 transition, which is similar to the U 1 to AFM ± transition peak of Fig. 1(a).

F. Alternating level shifts: Quantum fluctuations
The effects of the quantum treatment are similar to those found for uniform level shifts. The quantum system now always exhibits an AFM phase as the alternating level shifts explicitly break the translational symmetry. The ICR for the quantum model shows peaks around the mean-field phase transitions and we again find enhanced fluctuations in the excitation number around regions of bistability, but not in the photon number; Fig. 5. The enhanced fluctuations appear to agree much better with the mean-field contours, especially around the resonance of the checkerboard mode. Interestingly, there is a large decrease in fluctuations around the resonance of the checkerboard subradiant mode.

IV. DISCUSSION
In the limit of LLI, two-level atoms respond to light as linear classical oscillators [66]. Although atom-by-atom simulations of such systems, especially in large randomlydistributed ensembles with light-induced positions correlations between the atoms, can be demanding on numerical resources [62], the number of equations scales linearly with the atom number. Finding full quantum solutions in large systems, however, becomes quickly prohibitively challenging as the size of the density matrix in Eq. (1) scales exponentially with the atom number ∼ 2 2N . In this paper, we have approximated the quantum dynamics of a large array by subjecting it to periodic boundary conditions and varying parameters of only four atoms. This approach, however, provides a useful comparison with the corresponding mean-field dynamics of Eqs. (3) by unambiguously identifying quantum effects in the differences between the responses of the two cases.
While identifying light-established quantum correlations is interesting on its own right, this leads to practical implications as the number of equations in the meanfield dynamics scales linearly with the atom number. Determining the limits of validity of mean-field models can therefore provide a range of useful computational tools for the appropriate regimes. There is also a more philosophical point of view: As experiments with pristine quantum control of small atomic systems with genuine multimode dynamics are improving, the interface between quantum mechanics and classical physics, and the transition to classical physics due to decoherence or quantum stochastic nonlinear phenomena, is becoming ever more relevant in many-body systems. When a classical system exhibits the most dramatic consequences of nonlinearity, such as phase transitions or bistability, also the most recognisable differences between the quantum and classical theories arise. Instability in a classical phase transition represents exponentially growing deviations from the unstable solution to a new stable one, and bistability the simultaneous existence of two stable solutions. Quantum mechanics typically cannot favour either of the corresponding solutions, the dynamics are determined by the initial conditions and the evolution can also emerge as a superposition state, resulting in an enhanced fluctuations of measurement observables, as those identified in our study. ACKNOWLEDGMENTS J.R. and C.D.P. acknowledge financial support from the UK EPSRC (Grant Nos. EP/P026133/1, EP/M013294/1)

Appendix A MODEL FOR LIGHT-ATOM COUPLING
We express the electrodynamics in the length gauge, obtained by the Power-Zienau-Wooley transformation [67], such that E ± (r) = D ± F (r)/ 0 correspond to the positive and negative frequency components of the electric displacement in free space. The many-body quantum master Eq. (1) is then expressed in the rotating-wave approximation in terms of slowly varying field amplitudes and atomic variables, where E + e iωt → E + and σ − l e iωt →σ − l . The dipole-dipole interaction term is given by the real and imaginary part of the dipole radiation kernel, where γ jj = γ is the single-atom linewidth. The dipole radiation kernel acting on a dipole located at the origin yields the familiar dipole radiation expression [56] with r = |r|,r = r/r.
For the observables, such as transmitted light intensity or the photon count rate, the electric field product can be expanded in terms of incident and scattered fields to give whereÊ −Ê+ is the dyadic product with elements E α E * β , with α and β denoting the vector components. The first term in Eq. (14) gives the incident intensity, while the next three terms are the coherent scattered light, which remain even in the absence of quantum fluctuations. The last term, δÊ − s (r)δÊ + s (r ) = Ê − s (r)Ê + s (r ) − Ê − s (r) Ê + s (r ) , is the incoherent scattering, which is light scattered by disorder and quantum fluctuations. Because we consider atoms at fixed positions, the incoherent scattering is determined purely by quantum correlations. We measure the incoherent scattering from the ICR, obtained by substituting δÊ − s (r)δÊ + s (r ) into the photon count rate expression, Eq. (5), which gives The ICR that we will use under the mean-field description reads [20] As noted in the main text, this expression differs from the usual semiclassical description of the incoherent scattering [57] (which would vanish for fixed atomic positions) due to the inclusion of σ ee l terms. When calculating the photon count-rate, we integrate the field over a solid angle with NA = sin θ = 0.24, except when looking at photon fluctuations in Fig. 5, where we integrate over a closed surface.
Coherently transmitted light through a finite array can be approximated at a point (0, 0, ξ) from the centre of the array from [61,63,68,69] when λ ξ √ A, where A is the total area of the array and l is summed over all images of the atom j that are included due to the periodic boundary conditions. The total coherent field is thenÊ + (r) = 4 jÊ + j (r). We have found numerically this approximation works well.

Appendix B COLLECTIVE LOW LIGHT INTENSITY EIGENMODES
In the limit of LLI, the coupled-dipole model or the classical linear oscillator model becomes exact for the two-level atoms [37,66]. In the periodic lattice system, the LLI collective excitation eigenmodes are obtained by diagonalizing Eq. (12), and are the Bloch waves v (+) q (r l ) = A q cos(q · r l ), where A q = 2/N except for A q=0,(π/a,π/a) = 1/ √ N . The wavevectors q have components q x/y = 2πm x/y /N x/y a, where m x/y = 0, 1, .. to N x/y /2 or (N x/y − 1)/2 for an even or odd number of sites respectively, and N x/y is the number of sites along the x/y direction. The corresponding eigenvalues, δ q + iυ q , represent the collective line shifts (from the resonance of an isolated atom) δ q and linewidths υ q .
In our system there are four relevant LLI eigenmodes: the spatially uniform mode, the spatially nonuniform modes with striped phase variation alongx ±ŷ, and a checkerboard phase variation, q=(π/a,π/a) (r l ) = 1 2 cos π a , π a · r l . (22) To simulate an infinite system, we use periodic boundary conditions by adding repeat images of the system to the boundaries. We truncate to 101 images along thex and y direction, which gives an effective lattice size of 406 × 406. Numerically, the linewidths are given by υ uni = 25γ, υ + = 0.09γ, υ − = 0.08γ, and υ cb = 2 × 10 −4 γ. Due to our image truncation, the two striped modes υ ± are not degenerate and all 3 subradiant modes have a nonzero linewidth. The population of the eigenmodes [ Fig. 2(b)] is calculated using the occupation measure defined by [70] L α = N l v α (r l )ρ

Appendix C BIPARTITE ANSATZ
Most of the steady-state solutions to Eqs. (3) are obtainable using a bipartite lattice ansatz, where the 4 independent atoms are divided into 2 sites, labelled A and B. For spatially uniform level shifts, two divisions of the lattice can be chosen: sites A and B can lie along thex+ŷ andx−ŷ directions, respectively, forming a checkerboard variation, or both A and B sites lie along thex (ŷ) direction, alternating with respect to each other in theŷ (x) direction, forming a striped variation. For the nonuniform checkerboard level shift profile, the A and B sites must lie along thex+ŷ andx−ŷ directions, respectively, due to the spatial variation of the level shifts.
The steady-state solutions to Eqs. (3) under the bipar-tite ansatz obey where the subscript A/B denotes the lattice site group and we have made use of the symmetry under A ↔ B.