Quantum and nonlinear effects in light transmitted through planar atomic arrays

Understanding strong cooperative optical responses in dense and cold atomic ensembles is vital for fundamental science and emerging quantum technologies. Methodologies for characterizing light-induced quantum effects in such systems, however, are still lacking. Here we unambiguously identify significant quantum many-body effects, robust to position fluctuations and strong dipole–dipole interactions, in light scattered from planar atomic ensembles by comparing full quantum simulations with a semiclassical model neglecting quantum fluctuations. We find pronounced quantum effects at high atomic densities, light close to saturation intensity, and around subradiant resonances. Such conditions also maximize spin–spin correlations and entanglement between atoms, revealing the microscopic origin of light-induced quantum effects. In several regimes of interest, our approximate model reproduces light transmission remarkably well, permitting analysis of otherwise numerically inaccessible large ensembles, in which we observe many-body analogues of resonance power broadening, vacuum Rabi splitting, and significant suppression in cooperative reflection from atomic arrays. Many-body formalism for light transmission is developed to identify light-induced quantum effects and entanglement in atomic planar arrays and disks. Outside the regimes of pronounced quantum effects, potent semiclassical models reveal many-body analogues of power broadening and vacuum Rabi splitting.

L ight can mediate strong interactions between atoms, inducing strong position-dependent correlations, even in the limit of low light intensity, when the response (for the case of a simple level structure) is entirely classical. Such a correlated optical response can differ dramatically from that predicted by standard electrodynamics of continuous media, where resonantlight-induced dipole-dipole (DD) interactions between atoms are treated in an averaged sense 1,2 . Beyond the limit of low light intensity, an isolated atom can scatter light quantummechanically (such as in the incoherent Mollow spectrum 3 ), and quantum effects in the interactions of light with dilute atomic ensembles have been utilized in, e.g., quantum information protocols 4 . In strongly interacting dense systems the possible role of quantum and cooperative effects is less clear and has been the subject of long-standing debate 1,[5][6][7][8] . A particularly promising system to explore and utilize strong light-induced DD interactions is a regular planar array of scatterers, such as atoms. As shown both theoretically and experimentally, in the linear low-excitation limit these manifest a wealth of phenomena, e.g., subdiffraction features 9,10 , nontrivial topological phases 11,12 , transmission varying from complete reflection to full transparency [13][14][15][16][17][18][19] , narrow resonances, and subradiance [15][16][17][20][21][22][23][24][25][26][27] , as well as quantum technological applications 28,29 and other collective effects [30][31][32][33][34][35] .
We show that we can identify quantum effects in the light transmitted through planar arrays and uniform-density disks of cold and dense atomic ensembles. Many-body quantum correlations are induced by light-atom coupling, which, surprisingly, survive even strong many-body resonant DD interactions and atomic position fluctuations. Specifically, comparing the correlated optical response determined using the quantum master equation (QME) to simulations neglecting any quantum fluctuations between atomic levels in different atoms [referred to as the "semiclassical" equations (SCEs) 36 ], we systematically identify light-established quantum effects between atoms in the transmitted light as a function of atom confinement, density, and driving intensity. The effect of many-body quantum fluctuations on the scattering manifests most prominently at high densities when the light is close to saturation intensity, and especially significantly in the vicinity of subradiant resonances. We find that these conditions also produce maximal spin-spin correlations and entanglement of formation in the underlying atomic system, further confirming the role of many-body quantum correlations and entanglement in observing a difference in light transmission between QME and SCEs models. Incorporating the single-atom quantum description of light emission into the semiclassical scattering, we can typically use SCEs also for incoherent scattering to qualitatively reproduce the full quantum scattering even in the regimes where quantum effects in coherent scattering are most pronounced, and elsewhere also quantitatively. SCEs therefore allow us to analyze cooperative transmission of light through large atomic arrays and disks beyond the limit of low light intensity, without needing to solve the full strongly-interacting quantum dynamics. Doing so, we find collective phenomena due to DD interactions that are many-body analogs of power broadening and vacuum Rabi splitting of atomic resonances in cavities 37,38 , and demonstrate a significant effect of intensity on the transmission that may ultimately restrict the utilization of atomic arrays as highly reflective cooperative mirrors.

Results
An appealing feature of light scattering from cold atoms [39][40][41][42][43][44][45][46][47][48][49][50] is that light-mediated strong DD interactions can establish correlations between atoms at fluctuating positions, which are most simply described using atomic field operators for the ground and excited statesψ g;e ðrÞ. Hence, hP þ ðrÞi ¼ hψ where a, b, c, d ∈ {g, e}, representing the correlations in the optical response of an atom at r given the presence of a second atom at r 0 . These in turn depend on three-body correlations, etc., resulting in a hierarchy of equations of motion for correlation functions 51,52 . In a cold, dense ensemble ( Fig. 1) this hierarchy can significantly and nonperturbatively modify the scattering behavior even in the classical regime, invalidating attempts to truncate it 1 . This is a key ingredient in, e.g., Anderson localization of light, which has been a subject of considerable controversy and debate 53,54 .
Quantum master equation. A numerical device for solving this correlation function hierarchy is to treat the atoms as discrete point particles, meaning for a particular configuration of atomic positions {r 1 , …, r N } that two-body correlation functions take the formψ where ρ ðj;'Þ ad;bc denote correlation functions of the internal atomic energy levels only 36 . We then solve the internal atom dynamics at discrete positions, and the new correlation functions simply emerge from the N-body density matrix ρ ρ fr 1 ; ;r N g . This evolves according to QME: where the collective scattering is represented by the dispersive Ω jℓ and dissipative γ jℓ DD interactions, the single-atom half-width at half-maximum (HWHM) linewidth by γ jj = γ, and σ For simplicity, we consider two-level where E þ ¼ ðE À Þ Ã is the positive-frequency-component of the frequency ω = kc = 2πc/λ laser field, detuned from the atomic resonance frequency ω ge by Δ = ω − ω ge andσ ðjÞ ee ¼ e j i jj e h j. We take the polarization of the incident field to be parallel to the orientation of the atomic dipoles. Spatial correlations are numerically synthesized by ensemble-averaging over stochastic realizations of atomic positions sampled from the density distribution 36,55 . Solving Eq. (2) for large systems is numerically taxing, although few-atom ensembles already demonstrate manybody effects in their spectra 56 . Semiclassical model. In the limit of low light intensity, where the excited state population vanishes, the internal level correlations, such as those described by ρ ðj;'Þ ad;bc in Eq. (1), also vanish for twolevel atoms. The stochastic electrodynamics simulations are then formally exact 36,55 , reproducing the many-atom spatial correlations, which are identical to those occurring in the classical electrodynamics of coupled linear electric dipoles. Beyond the limit of low light intensity, the full dynamics of Eq. (2) can be greatly simplified by factorizing the internal atomic level correlation functions: Following the formalism of ref. 36 we then obtain coupled nonlinear equations: Note the relatively small number of equations 2N compared to the full quantum system size 2 N . This formalism has been applied to the modeling of pumping of atoms in dense clouds 50 , and has also been extended to cavity quantum electrodynamics (QED) 57 . Spatially correlated scattering between different atoms is accounted for in Eqs. (5) and (6) via Ω jℓ and γ jℓ (for Ω j≠ℓ = γ j≠ℓ = 0 they reduce to the independent-atom optical Bloch equations). In the limit of low light intensity the ensemble-averaged response of SCEs coincides with the exact classical electrodynamics; beyond this limit, the model incorporates nonlinear internal level dynamics of the atoms. However, because of the factorization in Eq. (4), they cannot account for many-body quantum entanglement between different atoms' internal levels. Finding situations in which the predictions of SCEs observably differ from the full QME solution therefore identifies light-induced quantum effects in the transmitted light. Conversely, regimes where quantum fluctuations are minimal allow for the simulation of much larger systems than are accessible with QME , and also test the validity of related approaches in other contexts, based, e.g., on mean-field approximations, intensity expansions, or truncations of the correlations [58][59][60][61][62] .
Light-established quantum effects in transmitted light. We begin by calculating the coherent (Fig. 2a-c) and incoherent ( Fig. 2d-f) forward transmission, T coh and T inc (see "Methods" section), through planar square arrays and thin disks of N = 4 atoms (Fig. 1), and the corresponding relative differences (Fig. 3a-c) between quantum and semiclassical results (for the effects of a larger 3 × 3 array, see Supplementary Note 1). The array could be realized, e.g., by an optical lattice 63 or dipole traps 64,65 . Unless otherwise stated, we consider lattice spacing a = 0.25λ and disk radius R = 0.28λ. Physically, we calculate the far-field light intensity in the same mode as the driving field E ± ðrÞ, integrated over the polar angle sin θ ≲ 0:24 (see "Methods" section and Supplementary Note 2). We account for the fluctuations in atomic positions due to finite trap confinement by ensemble-averaging over many stochastic realizations of position configurations (see "Methods" section) 10 . Light-induced positiondependent correlations between the atoms (see earlier "hierarchy of correlation functions") exist also classically, since the DD interactions depend on the precise interatomic separations. These classical many-body position correlations are therefore present in both the full quantum T QM coh and semiclassical T SC coh transmissions, but all quantum effects have been neglected in T SC coh . Hence, we unambiguously identify quantum effects in the coherent transmission by the difference T QM coh À T SC coh . Moreover, since the coherent scattering for a single atom is always purely classical, T QM coh À T SC coh cannot depend on the single-atom response, and therefore represents quantum correlations solely due to manybody effects.
To obtain the incoherently scattered light hδÊ 66 . This yields the incoherent transmission (see "Methods" section), from which we also isolate quantum behavior by T QM inc À T SC inc , in which case all the quantum effects have again been systematically neglected in T SC inc . Although the coherent scattering of a single atom is classical, this is not the case for the incoherent emission. We can improve the semiclassical incoherent model, without increasing the computational complexity, by adding the single-atom quantum description of incoherent light emission for all the atoms. In a single stochastic realization of atomic positions, the incoherent scattering contribution to intensity from independent quantummechanical atoms / P j A j ðhσ ðjÞ ee i À jhσ ðjÞ þ ij 2 Þ, where A j encapsulates the light propagation effects (see "Methods" section) 66 . Augmenting the semiclassical model with this single-atom quantum description integrated over the sample yields the incoherent transmission T SAQ inc . The quantum effects of the incoherent signal solely due to many-body processes are then encapsulated in T QM inc À T SAQ inc . The optical depth of the coherent transmission (Figs. 2a-c, 3a, c) corresponds physically to the degree to which extinction of the incident laser field by the averaged scattered field acts to reduce the transmission. The incoherent fluctuations in the scattered field (Figs. 2d-f, 3b) then counteract this reduction. In Fig. 2b, we identify many-body quantum fluctuations in the coherent transmission (T QM coh À T SC coh ) that increase with increasing DD interaction (Fig. 3a, c), reaching normalized residuals of over 10% at a = 0.25λ and I ≃ I sat (when the dipole amplitudes are greatest), where I ¼ 2ϵ 0 cjE þ ðr ¼ 0Þj 2 is the peak incident laser intensity and I sat = 4π 2 ℏγc/3λ 3 is the saturation intensity. Remarkably, even for a fully random disk quantum effects on the scattering do not wash out, but can produce residuals between the models of a few percent. On the other hand, there are also regimes where T SC coh accurately describes transmitted light. For example, in Fig. 3c the difference between T QM coh and T SC coh is less than 5% for a ≳ 0.4λ or I/I sat ≳ 16.
In Fig. 3c we observe two distinct peaks in the quantum manybody signatures of the coherent scattering for a = 0.25λ. The peak at ffiffiffiffiffiffiffiffiffiffi ffi I=I sat p % 0:3 results from a Fano interference between a narrow and an overlapping broad resonance (Fig. 3c inset), originating from an underlying highly subradiant (HWHM = 0.1γ) and superradiant (HWHM = 2.7γ) low light intensity eigenmode, respectively. The incident light couples to the phase-matched uniform superradiant eigenmode, but due to nonorthogonality of the non-Hermitian eigenmodes (Supplementary Note 3), the subradiant mode with rapid phase variation ("checkerboard" phase pattern, with dipoles at adjacent sites π out-of-phase) becomes populated at the narrow resonance (Fig. 2a). Strikingly, around this narrow resonance quantum effects constitute over 30% of the coherent scattering signal. We observe in Fig. 3c that, with increasing lattice spacing, the narrow Fano peak at ffiffiffiffiffiffiffiffiffiffi ffi I=I sat p % 0:3 disappears, as the corresponding eigenmode becomes less subradiant.
In Fig. 2e, f, on the other hand, we see the incoherent transmission is almost entirely dominated by quantum fluctuations when I ≳ I sat (T SC inc ! 0). However, once we incorporate the single-atom quantum description into the scattering and therefore transmission T SAQ inc , the difference becomes much smaller and the many-body quantum fluctuations are, as with the coherent scattering, maximal around I~I sat . Hence, using the improved model T SAQ inc , it is possible, even for incoherent scattering, to obtain excellent qualitative, and frequently quantitative agreement with the full quantum scattering (see Supplementary Note 4 for further demonstration of the different roles of quantum fluctuations in the incoherent scattering).
Spin correlations and entanglement. Up until now, we have identified many-body quantum effects via the transmitted light. These originate from the light-induced quantum correlations between internal levels of different atoms that do not satisfy the factorization assumption [given in Eq. (4)] of SCEs. In Fig. 4 we explicitly show these induced spin-spin correlations, and in Fig. 5a-c we show the many-body entanglement of formation (for a pair of atoms, using the formalism of ref. 67 -see Supplementary Note 5). The spin-spin correlations and entanglement exhibit behavior qualitatively similar to the quantum many-body correlations observed in the light scattering. As in Fig. 3c, both the correlations (Fig. 4) and entanglement (Fig. 5a-c) are maximal at I~I sat , where the intensity at which this peak occurs and the peak's amplitude decrease for increasing atomic spacing. The correlations and entanglement also both manifest linesplitting for I ≳ I sat , corresponding to when the atomic excited state population starts to increase. This is consistent with Fig. 2, where the observed quantum effects in transmitted light change from being maximal close to zero detuning when I~I sat , to being maximal off resonance when I ≫ I sat . The purity of the atomic state ( Fig. 5d) decreases to the limiting value 1/4 (Supplementary Note 5) for increasing I and is less sensitive to the array spacing than the entanglement E.
Large ensembles. When conditions are such that quantum effects on the light scattering are minimal, we can neglect quantum fluctuations, employing SCEs [Eqs. (5) and (6)] to analyze the coherent transmission through much larger ensembles, for which the full QME is inaccessible. In Fig. 6a-c we show how the transmission lineshapes of a 10 × 10 array differ significantly from the Lorentzian of independent atoms. For a single atom the linewidth is power broadened γ PB ðIÞ ¼ γ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 þ I=I sat p . For spacings a ≳ 0.3λ, the coherent optical depth for a 10 × 10 interacting array can also be fitted well to a single Lorentzian (ignoring the small structure present at low light intensity due to interfering eigenmodes, Fig. 6a), with a linewidth which also scales with intensity. If the linewidth at low light intensity is superradiant (e.g., a = 1.1λ) or subradiant (e.g., a = 0.8λ), it will also be correspondingly larger or smaller than γ PB for high intensities (Fig. 6c). For smaller spacings (e.g., a = 0.25λ), at high intensity the lineshape splits into a double resonance, producing a dip or "hole burning" (Fig. 6b). This dip is analogous to vacuum Rabi splitting 37,38 , where the interatomic DD coupling has now taken the role of the cavity coupling. In cavity QED, the mirrors create images of an atom inside the cavity, mimicking a periodic array, and the resonance doublet can be understood as a splitting of the excited state (Supplementary Note 6). While the dip in Fig. 6b only occurs at sufficiently high density, it can interestingly still exist even in the fully random ensemble. With increasing light intensity, the incoherent scattering lineshape tends to that of the independent atom (Supplementary Note 7). Fig. 3 Quantum many-body effects in optical transmission. a Relative differences, diffðSÞ maxfjS SCEs ðΔÞ À S QME ðΔÞj=S QME ðΔÞg, between quantum and semiclassical coherent optical depths (OD QM coh ¼ Àlog ðT QM coh Þ and OD SC coh , respectively) for a 2 × 2 array with lattice constant a = 0.25λ and normally distributed position fluctuations with varying rms width σ x,y = σ, σ z = 0.025λ (solid markers and thin solid lines). The open markers and thick solid lines are for a disk of randomly distributed atoms with radius R = 0.28λ and σ z = 0.025λ. The light intensity is ffiffiffiffiffiffiffiffiffi ffi I=I sat p ¼ 0:3 (circles), 2 (squares), and 10 (triangles). b Relative differences between quantum and semiclassical (augmented by a single-atom quantum description) incoherent transmission, T QM inc and T SAQ inc respectively, for the same parameters as in a. c Relative differences between quantum and semiclassical coherent optical depth as a function of light intensity, for fixed atomic positions and for different spacings. The inset shows OD QM coh (solid) and OD SC coh (dashed) for a fixed array with a = 0.25λ and ffiffiffiffiffiffiffiffiffi ffi I=I sat p ¼ 0:3. The incident beam has waist w 0 = 10λ, polarization ðx þŷÞ=  ad;bc respectively, for a 2 × 2 array with spacing a = 0.25λ and atom indices (j, ℓ), as illustrated in a. For fixed atomic positions, the correlations are of solely quantum origin, occurring between internal levels of different atoms. a The correlation lineshape for increasing intensity ffiffiffiffiffiffiffiffiffi ffi I=I sat p ¼ 1 (reddotted), 2 (blue-dashed), 5 (green-dot-dashed) exhibits an energy splitting. b The correlations peak around I~I sat , tending to zero at very low and high light intensities. The behavior is qualitatively similar for other spin/atom pairs. The incident beam has waist w 0 = 10λ and polarization ðx þŷÞ= ffiffiffi 2 p . A key feature of general subwavelength-spaced resonant emitter arrays is that they can exhibit perfect reflection 13,14 , which may typically be modeled using point-dipole scatterers [15][16][17]19 . Dipolar planar arrays can act as cooperative antennae 68 , with applications to quantum information processing 29 , making understanding nonlinear transmission essential. We calculate this for large arrays in Fig. 6d, and find that the reduction in the extinction as a function of light intensity is considerable-although less prominent with smaller spacings. This may ultimately restrict the application of atomic arrays as highly reflective cooperative mirrors to weak light intensities only.

Discussion
By comparing different descriptions of atom-light dynamics and light scattering, we have identified light-induced quantum effects, spin-spin correlations, and quantum entanglement in the light transmitted through planar arrays and disks which survive both position fluctuations and strong DD interactions. At narrow subradiant resonances, quantum fluctuations can be over 30%. Outside these resonances, provided we improve the model by incorporating the single-atom quantum description, SCEs typically still reproduce, also for incoherent scattering, the full quantum behavior at least qualitatively. This provides a methodology to calculate transmission of light through large arrays, consisting of hundreds of atoms, which can exhibit striking many-body phenomena (even without any quantum effects) reminiscent of single-atom power broadening and vacuum Rabi splitting. The existence of many-body quantum effects despite strong driving, high densities, and even with significant atomic position fluctuations is surprising. It suggests that optical quantum information processing in atomic ensembles 4 need not necessarily be restricted to dilute systems. Subradiant resonance narrowing has now been experimentally observed in the transmitted light through an optical lattice of atoms in a Mott-insulator state in the classical limit of low light intensity 69 . Several of our findings could also be verified in this setup by increasing the intensity of the incident light. The presence, even in uniform disks, of many-body effects attracting considerable interest 1,2,49,70 further relaxes the conditions necessary for their experimental observation.

Methods
Dynamics and correlation functions. We simulate the optical response of N-atom ensembles by stochastically sampling fixed positions {r 1 , …, r N } of stationary atoms, as the atomic center-of-mass dynamics are assumed negligible. In the full quantum dynamics, for each stochastic realization we solve the equations of motion for the Natom density matrix ρ fr 1 ; ;r N g ðtÞ with the atoms at fixed positions {r 1 , …, r N }, obeying QME [Eq. (2)]. In the Hamiltonian in Eq. (3) we use slowly-varying field amplitudes and atomic variables where the rapid rotation at the laser frequency ω has been factored out by substitutions E þ e iωt ! E þ ,σ ðjÞ À ðtÞ e iωt !σ ðjÞ À ðtÞ, etc. The collective coupling matrices Ω jℓ and γ jℓ , resulting, respectively, in collective resonance line shifts and linewidths in Eq. (2) (see Supplementary Note 3), are the real and imaginary parts of the dipole radiation kernel G(r): where is the electric field amplitude for an oscillating electric dipoled at the origin and n ¼ r=r. Note that we typically drop the contact interaction term 36 . Once ρ fr 1 ; ;r N g ðtÞ is known, the one-body ρ and so forth. Here ρ ðjÞ ge represents the correlations hψ y gψe i of Eq. (14), and corresponds to the matrix element 〈e|ρ 1 |g〉 of the one-body density matrix ρ 1 . Fig. 6 Semiclassically evaluated optical depth, resonance Rabi splitting, power broadening, and maximum extinction for 100 atoms. a, b From bottom to top: semiclassical coherent optical depth OD SC coh for a 10 × 10 square array with lattice constant a = 0.25λ and normally distributed position fluctuations (rms width σ x,y /a = 0, 0.1, 0.2, respectively); and 100 atoms randomly distributed within a radius R = 1.4λ (peak density 1.0k 3 ) thin disk. The standard sampling errors (shaded regions) are too small to see. For clarity each line is offset from that below it by 0.5. Fluctuations in z are characterized by a σ z = 0.025λ normal distribution, except when atoms are fixed (σ x,y = σ z = 0). The incident beam has waist w 0 = 10λ, polarization ðx þŷÞ= ffiffiffi 2 p , and intensity I/I sat = 10 −4 (a), and I/I sat = 100 (b). Black dashed lines show normalized lineshapes for a single atom with power-broadened γ PB ¼ γð1 þ I=I sat Þ 1=2 . c HWHM of OD SC coh for atoms at fixed positions with a = 0.8λ (squares), 1.1λ (circles), and for a single atom (black dashed line). The beam (w 0 = 30a) is nearly uniform over the atoms and the light is collected over a solid angle sin θ ≲ 0:04 (a = 0.8λ) and 0.03 (a = 1.1λ). d Peak extinction (1 À T SC coh ) with fixed positions: a = 0.8λ (squares), 1.1λ (circles), and 0.4λ (triangles). We set the beam waists to optimize extinction (3λ for a/λ = 0.8, 1.1 and 1.5λ for a/λ = 0.4) over a collection solid angle sin θ ≲ 0:37 using circular polarization Àðx þ iŷÞ= ffiffiffi 2 p .
In each stochastic realization, the N-atom configuration of positions {r 1 , …, r N } is obtained by sampling from a joint probability distribution P(r 1 , …, r N ), taken to be the initial distribution of stationary atoms. Ensemble-averaging over many such realizations then transforms the expectation values ρ and it is through solving the coupled dynamics between the light and atoms for each stochastic run and ensemble-averaging over many such realizations that we establish the light-induced spatial correlations between atoms 36,55 . For a single, isolated atom at the laser focus, the solutions to Eqs. (5), (6) (i.e., the optical Bloch equations) in the steady state are where I=I sat ¼ 2jd eg Á E þ ð0Þj 2 =ð_γÞ 2 . The effective linewidth in the denominators of Eqs. (16) and (17) is now power broadened, i.e., γ PB ¼ γ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 þ I=I sat p . We compare the full quantum solution of QME [Eq. (2)] with SCEs [Eqs. (5) and (6) In general for the atomic distribution before the light enters the sample we have P(r 1 , …, r N ) = |Ψ(r 1 , …, r N )| 2 , where Ψ(r 1 , …, r N ) denotes the N-body atomic wave function in position representation. For the initially uncorrelated atoms, each atom is sampled independently. We consider two different geometries: (i) atoms trapped in a two-dimensional array with precisely one atom per site; and (ii) a random, uniform distribution of atoms inside a thin, cylindrical disk of radius R and thickness Z. For the former case, we can sample the stochastic position of an atom in each site 10 , obtaining P(r 1 , …, r N ) = ϱ 1 (r 1 )…ϱ N (r N ), where the density distribution of the jth array site, centered at r j = R j , is approximated by a Gaussian: where the standard deviations σ x , σ y , σ z quantify the spatial confinement of the trapped atoms in all three directions.
Scattered light. The total electric field operatorÊ ± ðrÞ ¼ E ± ðrÞ þÊ ± d ðrÞ is the sum of the laser field and the fields scattered from all atoms To analyze the different contributions in the scattered light, we write it aŝ E We then obtain, hereÊ ÀÊþ is a dyadic product with elementsÊ À αÊ þ β , with α, β ∈ {1, 2, 3} cycling over the different polarization components, where the intensity is proportional to its diagonal elements. The first term on the right hand side of Eq. (21) yields the incident light intensity, the second, third, and fourth terms produce the coherent scattering, and the final term produces incoherent scattering dependent on fluctuations. Rearranging Eq. (21) to solve for the incoherent scattering gives: δÊ which describes correlations in the scattered light.
Consider first a single atom at the origin R ¼ 0. According to Eq. (20), the coherently scattered light consists of expectation valueŝ and there is no difference between the quantum and semiclassical coherent scattering. Hence, any difference between the quantum and semiclassical coherent scattering for a many-atom ensemble is due solely to many-body quantum effects. The incoherent contribution in Eq. (22) is more subtle, aŝ means the incoherently scattered light from a single atom yields where we have usedσ þσÀ ¼σ ee . In the semiclassical approximation, where the quantum fluctuations are ignored, one then replacesσ þ by hσ þ i in Eq. (24) 66 , such thatÊ and the incoherently scattered light intensity in Eq. (25) vanishes. Unlike the coherent scattering, the incoherent scattering for a single atom therefore differs depending on whether we treat it in a quantum or semiclassical manner. Generalizing to the many-atom case, Eq. (24) now becomeŝ where, as in Eq. (20), ½Gðr À RÞ Ã acts onP À ðRÞ and likewise Gðr À R 0 Þ on P þ ðR 0 Þ. When calculating the full quantum solution, the correlation functions are evaluated using the solution to QME [Eq. (2)] and by ensemble-averaging over many realizations of atomic positions. However, we can also introduce the manybody version of the single-atom semiclassical approximation [Eq. (26) Deriving the semiclassical scattered light in Eq. (29) corresponds to a systematic way of neglecting all quantum fluctuations when the atomic response is first calculated from SCEs [Eqs. (5) and (6)]. Hence, comparing the scattered light of Eq. (29) with the equivalent full quantum solution of Eq. (27) provides a signature for quantum effects in the collective atomic response. Alternatively, if our goal is to determine a computationally efficient and accurate approximation to the full quantum solution, we can instead try to improve the semiclassical approximation. A simple way to achieve this without increasing computational demands is to include the single-atom quantum description to incoherent scattering [Eq. (25)] integrated over the extent of the sample, which is sufficient in a number of cases to capture the leading quantum contributions.
We begin this procedure by placing the atomic operators in Eq. (27) in the normal order. This yields for the expectation term on the right hand side of This improved description includes both the semiclassical contribution and the single-body quantum fluctuations, meaning any difference in the incoherent scattering between this improved model and the full quantum model is solely due to many-body quantum effects.
Transmission. Coherently and incoherently transmitted light is calculated through a disk of cross-sectional area S perpendicular to the optical axis a distance f = 500λ downstream of the atoms. We consider light transmitted in the same spatial mode as the driving field, motivated by a typical experimental scheme of collecting transmitted light into a single-mode optical fiber, although, for simplicity, we ignore any explicit refocussing or fiber coupling. The transmitted light therefore has the form Note that because of the double integral over S and S 0 the expectation term is now a function of r and r 0 , although substituting r 0 into the preceding equations does not affect our discussion of coherent and incoherent scattering.
To calculate the coherent transmission T coh (plotted as optical depth OD coh Àlog T coh ), we substitute the first four terms on the right hand side of Eqs. (21) into (33). This gives quantum T QM coh or semiclassical T SC coh coherent transmission, depending on whether we use the solutions to QME or SCEs . To calculate the incoherent contribution to the transmission, we replace the two-field expectation in Eqs. (33) with (22). Evaluating Eq. (22) using Eqs. (27) and (30), along with the solutions to QME , results in the quantum incoherent transmission T QM inc . Using instead the solutions to SCEs and either Eqs. (29) or (31), respectively, produces the semiclassical incoherent transmission T SC inc , or the improved model for incoherent transmission T SAQ inc where the independent-atom quantum description is added to the semiclassical model.

Data availability
The data presented in this study are available in the Durham University repository with the identifier https://doi.org/10.15128/r1zc77sq127. The codes used to generate these data are available in the Zenodo repository with the identifier https://doi.org/10.5281/ ZENODO.3924698.