Signatures of optical phase transitions in superradiant and subradiant atomic arrays

Resonant light interacting with matter supports different phases of a polarisable medium, and optical bistability where two phases coexist. Such phases have previously been actively studied in cavities. Here, we identify signatures of optical phase transitions and optical bistability mapped onto scattered light in free-space planar arrays of cold atoms. Methods on how to explore such systems in superradiant and extreme subradiant states 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 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. Non-equilibrium phase transitions induced by light have been of great interest in cavities where feedback is provided by cavity mirrors. Here, the fingerprints of phase transitions and intrinsic optical bistability, supported by dipole-dipole interactions alone, are identified in light scattered from planar atomic arrays.

R esonant 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 . In arrays of closely-spaced cold atoms, the strong light-mediated dipole-dipole interactions arise naturally, as atoms do not absorb light, their resonances are well defined, and the atoms can respond to light quantummechanically. Atomic arrays have been proposed as constituents of metamaterials 9 , for quantum information processing [10][11][12] , atomic clocks [13][14][15] , emission of nonclassical light [16][17][18] and entanglement [19][20][21] , and a way to realise topological phases 22,23 . In most experiments so far, the efforts to observe collective optical responses of cold atoms, in both random atomic ensembles [24][25][26][27][28][29][30][31][32][33][34] and in arrays 5 and chains 35 , 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 36 . 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 37 in atomic systems. Optical bistability and phase transitions have been actively studied in systems without the spatial correlations and structure of the sample [38][39][40][41][42][43][44][45] , e.g., in cavities where the feedback mechanism is provided by the cavity mirrors [46][47][48] . Intrinsic bistability is a process where phase transitions are generated by the self-interactions of the sample, and despite having been observed in highly-excited Rydberg atoms in the microwave regime 49 , 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 short-and long-range interactions [50][51][52] .
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 mean-field 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-fieldtheoretical 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 can be seen by peaks and angular dependence of the incoherent photon count rate, and most clearly in fluctuations 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.

Results
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. 1a. We take the polarisation and the direction of the atomic dipoles to beê ¼ Àðx þŷÞ= ffiffi ffi 2 p 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 N = 4 atoms in a square array. However, we also discuss Fig. 1 A 2D array of atoms illuminated by incident light and resulting stable phases. a The array has a lattice spacing a = 0.1λ and is illuminated by an incident field, E, with only the central four atoms (green) independent, while the remaining array is obtained by periodic boundary conditions. Calculations in a system with nine (blue and green) instead of four independent atoms result in qualitatively similar phases. Phases for b uniform and c alternating profiles of level shifts as a function of incident intensity, I/I sat , and laser frequency detuning from the atomic resonance, Δ/γ. Spatially uniform (U), antiferromagnetic (AFM), oscillatory (OSC) and phases that are not uniform or AFM in nature (BD) emerge.
briefly the system with a N = 9 atom variation to study finite-size effects from the N = 4 system. 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 53 where the square brackets represent a commutator andσ þ j ¼ e j i jj g h j ¼ ðσ À j Þ y the raising operator, where e j i j and g j i 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 Methods for details). The Hamiltonian is given bŷ whereσ ee l ¼σ þ lσ À l , Δ l ¼ ω À ω ðlÞ 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 cjE 0 j 2 in units of the saturation intensity, I sat = ℏc4π 2 γ/3λ 3 , or the Rabi frequency, 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 factorisation of internal level correlations, hσ α where ρ ðlÞ ge ¼ Trfσ À lρ ðtÞg and ρ ðlÞ ee ¼ Trfσ ee lρ ðtÞg. We will use Eqs. (3) and (4) to determine the long-time phases and optical bistability that can occur in the system.
where Gðr À r l Þ is the dipole radiation kernel (Eq. (16) in Methods). We will compare the optical responses obtained from the full quantum dynamics of Eq. (1) with those calculated from the mean-field Eqs. (3) and (4). We consider coherently transmitted light in the forward direction, T coh ¼ jê Á hÊ À ðrÞij 2 =ê Á E À ðrÞ j j 2 , expressed in terms of the optical depth OD ¼ Àln ðT coh Þ. We also calculate the rate of scattered photons  (18) in Methods). The ICR that we will use under the mean-field description (Eq. (19) in Methods) is different from the usual semiclassical approximation for the incoherent scattering. The atom-light dynamics is solved from the mean-field Eqs. (3) and (4), but the single-atom quantum description of emitted light hσ ee l i is now included 20 for the scattered light, and the ICR no longer vanishes for atoms at fixed positions 54 .
Collective low light intensity eigenmodes. In the limit of LLI, ρ ðlÞ ee ¼ 0, and the mean-field Eqs. (3) and (4) 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 , where q denotes the wavevector of the LLI eigenmodes. Collective modes with υ q > γ (υ q < γ) are termed superradiant (subradiant). The two-level transition resonance wavelength defines the light cone, ωa/c = 0.2π, where modes with |q|a > 0.2π are completely dark in an infinite lattice.
We focus on two LLI eigenmodes (see Methods): the spatially uniform superradiant mode v un (r l ) (υ un = 25γ) and a subradiant mode with a checkerboard-patterned phase variation with every atom oscillating π out-of-phase 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. 2a. The checkerboard eigenmode v cb is particularly interesting as it is outside the light cone for small lattice spacings. In the following, we propose a protocol to excite v cb and, in principle, prepare coherent strongly subradiant excitations that exist outside the light cone. This can be achieved by breaking and restoring the symmetry with ac Stark shifts 55 of lasers (or microwaves) that form a checkerboard pattern of atomic level shifts from a standing-wave, cos 2 ½5kðx þ yÞ= ffiffi ffi 2 p , with the intensity varying along the lattice diagonalx þŷ and the intensity maxima separated by ffiffi ffi 2 p a. 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 56 , which can have a resonance wavelength of λ ≃ 2.6 μm and spacing of 206.4 nm, resulting in the effective lattice spacing a ≃ 0.08λ. Figure 2a, b shows how checkerboard-patterned alternating level shifts of 2γ 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. 2b). The ac Stark shifts break the symmetry of the lattice, which allows the checkerboard mode to couple to the incident field. Upon removing the level shifts, the symmetry is restored and the subradiant eigenmode completely decouples from the incident field again, but there is now an excitation stored in the subradiant mode outside the light cone.
Analytic results for optical bistability. Classifying the steady states of the mean-field solutions of Eqs. (3) and (4) 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 ρ ðlÞ ge ¼ ρ ge , ρ ðlÞ ee ¼ ρ ee , and Δ l = Δ into Eqs. (3) and (4). We then obtain the stationary states where we have defined which, withΩ ¼ P j≠l Ω jl andγ ¼ P 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. Equations (7) and (8) are equivalent to the familiar solutions of the independent-atom optical Bloch equations (Eqs. (30) and (31) in Methods), but with the Rabi frequency R replaced by R eff . As ρ ge appears on the both sides of Eq. (7) via R eff , we generally have multiple solutions. For two different coexisting stable solutions, we have optical bistability.
We can eliminate from Eqs. (7) and (9) the atomic variables and obtain an equation for the incident light field R ¼ RðR eff Þ. The bistability threshold is then found when djRj 2 =djR eff j 2 ¼ 0. This gives a cubic polynomial for jR eff j 2 (see Eq. (33) in Methods) 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 ¼ jRj 2 =ðΔ 2 þ γ 2 Þ, and the cooperativity parameter, The two solutions represent very different responses to the incident light. The first solution (10), termed the cooperative solution (in an analogy with the terminology of optical bistability in cavities 40 ), 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 solution (11), termed the single-atom solution, exists when p unsat > 2ðjCj þ Re½CÞ, and arises when the atoms react to the incident light almost independently, with R % R eff when jRj ! 1. 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 satisfỹ Ω 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 Supplementary Note 1). 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 57 (see also ref. 58 ) 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 [38][39][40][41][42]44,45 can now be most easily illustrated, and our adapted terminology motivated, at the specific value of Δ=γ 1 Ω=γ for which C ¼γ=2γ in Eq. (12) is real. The expression for the incident light field (Eq. (32) in Methods) 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 40 .
In cavity quantum electrodynamics (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~1equivalent 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 36,59,60 . 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. Figure 3a shows the upper and lower intensity thresholds for bistability as a function of lattice spacing, along with the bistability region and solutions in Fig. 3b, c. 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, 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.
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) and (4) for ka = 0.2π. Here the coupling of atoms is to the full free-space electromagnetic spectrum, but different optically-induced phases can also occur in systems with a periodic lattice where the light-atom coupling is just to a single mode of a cavity 61,62 , and phase transitions could additionally be induced, e.g., by feedback 63 .
In general, we find the system can exhibit spatially uniform phases, antiferromagnetic (AFM) phases and persistent oscillations (OSC), as well as different optical bistabilities. The detailed phase diagrams in Fig. 1b, c are calculated for a square array of N = 4 atoms with periodic boundary conditions. However, analogous behaviour is anticipated to emerge for different independent atom numbers, with additional phases due to the presence of more LLI eigenmodes. We have explicitly simulated a square array of N = 9 atoms with periodic boundary conditions and found a qualitatively similar phase diagram, with spatially uniform phases, OSC phases and regions of bistability, but instead of an AFM phase, we have found a spin density wave (SDW) phase with three-site periodicity, due to the obvious constraints of the lattice length. Transitions between different phases can result in small dips in the OD, and lead to large peaks in the ICR. Phase bistabilities are identified in large jumps in the OD and ICR, as well as hysteresis upon varying the laser frequency. For larger independent atom numbers, the AFM phases will re-emerge and coexist with the SDW phase, with additional SDW phases that have different periodicities as larger scale fluctuations are allowed in the system. We have found the region of U 1 /U 2 bistability to remain the same without any additional fluctuation-induced phase instabilities. Therefore, our analysis of two independent atom numbers indicates our results presented here are generic features of large lattices, i.e., peaks and dips in the ICR, OD and excited state population around meanfield phase transitions and bistabilities.
We first consider the case where the atomic level shifts are all equal; Fig. 1b. Beyond the obvious phases representing uniform low and high excitation numbers, labelled U 1 and U 2 (the q = 0 case of Eq. (27) in Methods), respectively, we interestingly also find stable phases with spontaneously broken translational symmetries and regions of two coexisting stable phases. Different phases can be distinguished from one other via the staggered level population order parameter, m ¼ P j e iQÁr j ð2ρ ðjÞ ee À 1Þ, where each phase has a specific value of Q which will maximise m and identify the phase, e.g., Q = 0 for the spatially uniform phase. The coherent and incoherent optical responses, OD and ICR, from the mean-field analysis are shown in Fig. 4a-d. While the uniform phases U 1 and U 2 vary smoothly into one another (white regions of Fig. 1b), 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. (32) in Methods). However, differences occur due to one of the uniform phases becoming unstable at positive detunings. The resonances of U 1,2 for I/I sat = 100 in Fig. 4a-d are both broad and correspond to the superradiant LLI uniform excitation eigenmode (Eq. (23) in Methods) at low intensities, appearing at the detuning Δ = −25γ. For the U 1 /U 2 bistability at I/I sat = 200, we find hysteretic behaviour upon sweeping the resonance from either red-or bluedetuned 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. 1b). These AFM phases appear at detunings resonant with the LLI excitation eigenmodes u ±,cb (Eqs. (24) and (25) in Methods), 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γ. They can be distinguished by the staggered level population order parameter, m, with Q = (π/a, π/a) for the AFM cb phase, and Q = (π/a, 0), (0, π/a) for the AFM ± phase, which reflect the symmetry of the corresponding LLI eigenmodes. 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. 4c 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. In general, we find the peaks for the AFM phases are always far narrower than those for uniform phases, which distinguishes the uniform and AFM phases in the ICR. No clear signature of this transition can be seen in the OD. We find that the AFM cb phase is bistable with the U 1 phase (yellow regions of Fig. 1b). Switching off the incident drive, the uniform phase decays superradiantly, and the AFM cb 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. This also explains why AFM phases result in narrower peaks than uniform phases due to their underlying subradiant nature. The hysteresis associated with the bistability could be utilised as a possible method for preparing subradiant excitations, when a steady-state superradiant mode is transformed into a subradiant one by a laser frequency sweep.
Both AFM phases also can become unstable (via Hopf bifurcations) resulting in an OSC phase. 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, and this allows the OSC phase to be distinguished from the AFM phases as the signal from the OSC phase will vary over a typical timescale of the order τ~5-10/γ. 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.
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. 4a-f 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 is incorporated in the analysis of the mean-field dynamics. However, the ICR still shows a resonance near the detunings where mean-field AFM transitions occur. Furthermore, this resonance depends on the collection angle of the photons (see Supplementary Note 2), while the broad resonance in the ICR where there are uniform phases in the mean-field shows no angular dependence.
It is generally known from past bistability studies [43][44][45]48,52,64 that mean-field bistabilities coincide with enhanced quantum fluctuations, which can be understood as tunnelling between the two mean-field solutions. The corresponding quantum distribution is then bimodal. The calculated incoherently scattered photon number fluctuations IoD n ¼ ðhn 2 i À hni 2 Þ=hni, wheren is the operator form of Eq. (6) with all the light collected over a closed surface, however, shows no signatures of enhanced fluctuations (Fig. 5a) but we find that the fluctuations of the excited level population (Fig. 5b), 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 mean-field solutions, and could be detected by resonantly transferring excited atoms to another level. This behaviour is also present for N = 9 independent atoms (see Supplementary Note 3) for U 1 /U 2 and SDW/U 1 bistabilities, and indicates that enhanced fluctuations around regions of bistability should be a general feature for larger independent atom numbers.
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 (with removed level shifts) represents subradiant checkerboard eigenmode (Eq. (25) in Methods) existing outside the light cone. Here we are interested in a continuous driving of the system, such that we will not restore the symmetry by removing the level shifts after exciting the mode. A weak level shift of 2γ is maintained to couple the eigenmodes and even the LLI eigenmode outside the light cone is radiating, instead of completely trapping the excitation. We now analyse how alternating level shifts influence 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. 1b now transform to checkerboard AFM cb phases in Fig. 1c, but can still be distinguished as having low and high excitations, labelled AFM 1 and AFM 2 . The coherent and incoherent optical responses, OD and ICR, from the mean-field analysis are shown in Fig. 6a-d. 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 non-uniform level shifts, the intensities at which non-trivial phases emerge are lower due to the checkerboard subradiant mode (Eq. (25) in Methods) 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 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. 2a. 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 the AFM 2,1 phase, respectively. For the BD 1 (BD 2 ) phase, the atoms alongx þŷ are in-phase (out-of-phase) with respect to each other, while the atoms alongx Àŷ are out-ofphase (in-phase) with respect to each other, and have different dynamics to the atoms atoms alongx þŷ. 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 1 phase remains. BD 1,2 phases occur because of the striped subradiant modes (Eq. (24) in Methods), 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. 1c for the AFM 1 to BD 1 transition, which is similar to the U 1 to AFM ± transition peak of Fig. 1b. 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 cb 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.

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 randomly-distributed ensembles with light-induced positions correlations between the atoms, can be demanding on numerical resources 59 , 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) and (4) by unambiguously identifying quantum effects in the differences between the responses of the two cases. Simulations of nine atoms give similar results to the four-atom case, but with the emergence of an SDW phase with higher periodicity due to an increasing number of LLI modes in the system. The optical signatures of the SDW phase is similar to the AFM phases, with narrow peaks in the ICR and dips in the OD, so our analysis of the four-atom system provides key insight into the dynamics and signatures of large systems.
While identifying light-established quantum correlations is interesting on its own right, this leads to practical implications as the number of equations in the mean-field 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. Instead 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.

Methods
Model for light-atom coupling. We express the electrodynamics in the length gauge, obtained by the Power-Zienau-Woolley 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 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. (17) 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, hδÊ 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 54 (which would vanish for fixed atomic positions) due to the inclusion of hσ ee l i terms. When calculating the photon countrate, we integrate the field over a solid angle with NA ¼ sin θ ¼ 0:24, except when looking at photon fluctuations in Fig. 5a, c, 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 58,60,68,69 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Þ ¼ P 4 jÊ þ j ðrÞ. We have found numerically this approximation works well.
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 36,66 . In the periodic lattice system, the LLI collective excitation eigenmodes are obtained by diagonalising Eq. (15), and are the Bloch waves v ðþÞ q ðr l Þ ¼ A q cosðq Á r l Þ; ð21Þ v ðÀÞ q ðr l Þ ¼ A q sinðq Á r l Þ; ð22Þ . 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, v un ðr l Þ v ðþÞ q¼0 ðr l Þ ¼ the spatially nonuniform modes with striped phase variation alongx ±ŷ, v ± ðr l Þ v ðÀÞ q¼ðπ=2a; ± π=2aÞ ðr l Þ ¼ and a checkerboard phase variation, v cb ðr l Þ v ðþÞ q¼ðπ=a;π=aÞ ðr l Þ ¼ 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 thê x andŷ direction, which gives an effective lattice size of 406 × 406. Numerically, the linewidths are given by υ un = 25γ, υ + = 0.09γ, υ − = 0.08γ, and υ cb = 2 × 10 −4 γ. Due to our image truncation, the two striped modes υ ± are not degenerate and all three subradiant modes have a nonzero linewidth. The population of the eigenmodes (Fig. 2b) is calculated using the occupation measure defined by 70 Analytic mean-field solutions. For a general form R l ¼ Re iqÁr l , where q is the wavevector of the drive, a solution to Eqs. (3) and (4) is given by ρ ðlÞ ge ¼ ρ ge e iqÁr l , with are the Fourier transforms of the real and imaginary parts of the dipole kernel, Eq. (16), respectively (excluding the self-interaction j = l). The number of excitations ρ ee obeys the following cubic equatioñ γðqÞ 2 þΩðqÞ 2 Â Ã ð2ρ ee À 1Þ 3 þγðqÞ 2 þΩðqÞ 2 À 2ΔΩðqÞ À 2γγðqÞ Â Ã ð2ρ ee À 1Þ 2 þ Δ 2 þ γ 2 þ 2jRj 2 À 2ΔΩðqÞ À 2γγðqÞ Â Ã ð2ρ ee À 1Þ þ Δ 2 þ γ 2 ð Þ¼0: ð29Þ We use the solutions to Eqs. (27) and (29) to describe the spatially uniform solutions and their bistability in the phase diagram with uniform level shifts. In the absence of incident light, Eq. (29) admits only one real solution of ρ ee = 0. For large intensities, the interaction terms in Eq. (29) become negligible and the coherence and number of excitations become identical to the noninteracting solutions to the optical Bloch equations with familiar power-broadened linewidths where we have used 2jRj 2 =γ 2 ¼ I=I sat . Other steady state solutions to Eqs. (3) and (4) can be found numerically that obey a bipartite ansatz, with more details in Supplementary Note 4.
Bistability threshold. The regimes of bistability between the spatially uniform phases can be predicted analytically for certain detunings. Without the loss of generality in the derivation, we can set in the following q = 0 for the drive. Substituting ρ ðlÞ ge ¼ ρ ge and ρ ðlÞ ee ¼ ρ ee into Eqs. (3) and (4) gives the steady-state solutions (7) and (8) that are equivalent to the solutions of the optical Bloch Eqs. (30) and (31), but with the Rabi frequency R replaced by the total external electric field on an atom (the incident field plus the scattered light from all the other atoms) in the array, R eff (Eq. (9)). By substituting Eq. (7) into Eq. (9), we obtain where C is given by Eq. (12) forΩ Ωð0Þ andγ γð0Þ. Equation (32) is also valid for a general plane-wave Rabi drive (q ≠ 0) if 0 → q and ρ ge ! ρ ge e iqÁr j in Eq. (9). Analytic solutions to Eq. (32) can be found by either ignoring the single R eff term, or the Δ 2 + γ 2 term in the denominator, which gives the cooperative (Eq. (10)) and single-atom (Eq. (11)) solutions in the main text, respectively. It is worth noting that for real C, Eq. (32) has similar form as equations determining bistability in cavities 40 , as discussed in the main text.

Code availability
The code that supports the findings of this study is available upon reasonable request.