Quantum and Nonlinear Effects in Light Transmitted through Planar Atomic Arrays

We identify significant quantum many-body effects, robust to position fluctuations and strong dipole--dipole interactions, in the forward light scattering from planar arrays and uniform-density disks of cold atoms, by comparing stochastic electrodynamics simulations of a quantum master equation and of a semiclassical model that neglects quantum fluctuations. Quantum effects are pronounced at high atomic densities with the light close to saturation intensity, and especially at subradiant resonances. We find an enhanced semiclassical model with a single-atom quantum description provides good qualitative, and frequently quantitative agreement with the full quantum solution. We use the semiclassical simulations for large ensembles that would otherwise be numerically inaccessible, and observe collective many-body analogues of resonance power broadening and vacuum Rabi splitting, as well as significant suppression in cooperative reflection from atom arrays.

We identify significant quantum many-body effects, robust to position fluctuations and strong dipole-dipole interactions, in the forward light scattering from planar arrays and uniform-density disks of cold atoms, by comparing stochastic electrodynamics simulations of a quantum master equation and of a semiclassical model that neglects quantum fluctuations. Quantum effects are pronounced at high atomic densities with the light close to saturation intensity, and especially at subradiant resonances. We find an enhanced semiclassical model with a single-atom quantum description provides good qualitative, and frequently quantitative agreement with the full quantum solution. We use the semiclassical simulations for large ensembles that would otherwise be numerically inaccessible, and observe collective many-body analogues of resonance power broadening and vacuum Rabi splitting, as well as significant suppression in cooperative reflection from atom arrays.
Light 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 resonant-light-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 quantum-mechanically, and quantum effects in the interactions of light with dilute atomic ensembles have been utilized in, e.g., quantum information protocols [3]. In strongly interacting dense systems the possible role of quantum and cooperative effects is less clear and has been the subject of long-standing debates [1, 4-7]. One particularly promising system to explore and utilize strong light-induced DD interactions is a regular planar array of scatterers such as atoms. In the linear low-excitation limit these manifest, as shown both theoretically and experimentally, a wealth of phenomena, e.g., subdiffraction features [8-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][36].
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) [37]), 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. 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 such regimes, 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 analogues of power broadening and vacuum Rabi splitting of atomic resonances in cavities [38,39], and demonstrate a significant effect of intensity on the transmission that may ultimately restrict the utilization of atomic arrays as highly reflective cooperative mirrors.
An appealing feature of light scattering from cold atoms [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, P + (r) = ψ † g (r) d geψe (r) denotes the light-induced atomic polarization, where d ge = d * eg is the dipole matrix element, and the populations are ψ † g (r)ψ g (r) and ψ † e (r)ψ e (r) . Because of the DD interactions, the polarization and populations also depend on two-body correlations ψ † a (r)ψ † b (r )ψ c (r )ψ d (r) , 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 . These in turn depend on three-body correlations, etc., resulting in a hierarchy of correlation function equations of motion [51,52]. In a cold, dense ensemble this hierarchy can significantly and nonperturbatively modify the scattering behavior, 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].
A numerical trick for solving this correlation function hier-archy is to treat the atoms as discrete point particles, meaning that for a particular configuration of atomic positions {r 1 , . . . , r N } the two-body correlation functions take the form where ρ ( j, ) ad;bc denote correlation functions of the internal atomic energy levels only [37]. We then solve the internal atom dynamics at discrete positions, and the new correlation functions simply emerge from the N-body density matrix ρ ≡ ρ {r 1 ,...,r N } . This evolves according to QMĖ where the collective scattering is represented by the dispersive Ω j and dissipative γ j DD interactions, the single-atom halfwidth at half-maximum (HWHM) linewidth by γ j j = γ, and For simplicity, we consider two-level atoms and the Hamiltonian where E + = (E − ) * is the positive-frequency-component of the frequency ω laser field, detuned from the atomic resonance frequency ω ge by ∆ = ω − ω ge andσ ( j) ee = |e j j e|. Spatial correlations are numerically synthesized by ensemble-averaging over stochastic realizations of atomic positions sampled from the density distribution [37,55]. Solving Eq. (2) for large systems is numerically taxing, although few-atom ensembles already demonstrate many-body effects in their spectra [56].
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 two-level atoms. The stochastic electrodynamics simulations are then formally exact [37,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. [37] we then obtain coupled nonlinear equations Note the relatively small number of equations 2N compared to the size of the full quantum system 2 N . This formalism has been applied to the modeling of pumping of atoms in dense clouds [57], and has also been extended to cavity QED [58]. 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 ensembleaveraged 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 [59][60][61][62][63].
We begin by calculating the coherent and incoherent forward transmission, T coh and T inc [64], through planar square arrays and thin disks of N = 4 atoms (Figs. 1 and 2). The array could be realized, e.g., by an optical lattice [65] or dipole traps [66,67]. Unless otherwise stated, we consider lattice spacing a = 0.25λ and disk radius R = 0.28λ. Physically, we calculate the far-field light intensity I in the same mode as the driving field E ± , integrated over the polar angle sin θ 0.24 [64]. We account for the fluctuations in atomic positions due to finite trap confinement by ensemble-averaging over many stochastic realizations of position configurations [10, 64]. We unambiguously identify quantum effects in coherent transmission from the difference between the full quantum and semiclassical transmissions T QM coh − T SC coh . Since the coherent scattering of a single atom is classical, this difference is due solely to many-body quantum correlations in the atomic response.
To obtain the incoherently scattered light δÊ − d (r) δÊ + d (r) , we write the scattered light field asÊ + d = Ê + d + δÊ + d , where δÊ + d denotes the fluctuations [68]. This yields incoherent transmission [64] for which quantum behavior also is isolated by T QM inc − T SC inc . We can improve the semiclassical incoherent model, however, without increasing the computational 1 complexity, by adding the single-atom quantum description of incoherent light emission for all the atoms. In a single realization of stochastic atomic positions, the incoherent scattering contribution to intensity from independent quantummechanical atoms ∝ j A j ( σ ( j) ee − | σ ( j) + | 2 ), where A j encapsulates the light propagation effects [64,68]. Augmenting the semiclassical model with this single-atom quantum description integrated over the sample yields the incoherent transmission T SAQ inc . The many-body quantum effects of the incoherent signal are then encapsulated in T QM inc − T SAQ inc . In Fig. 1(b), we identify many-body quantum fluctuations in the coherent transmission (T QM coh − T SC coh ) that increase with increasing DD interaction [ Fig. 2(a,c)], reaching normalized residuals of over 10% at a = 0.25λ and I I sat (when the dipole amplitudes are greatest). Strikingly, quantum effects constitute over 30% of the signal in the vicinity of the narrow subradiant resonant shown in the inset of Fig. 2(c). This may be due to the enhanced dipole magnitude of subradiance [16], the antisymmetry of the collective dipolar eigenvector [see diagram in Fig. 1(a)], or the rapidly varying Fano interference. It is remarkable that even for a fully random disk quantum effects on the scattering are not washed out but can produce residuals between the models of a few percent.
Conversely, as in Fig. 1(e,f), once I I sat the incoherent transmission is almost entirely dominated by quantum fluctuations (T SC inc → 0 [69]). 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. For example, in Fig. 2(c) the difference between T QM coh and T SC coh is less than 5% for a 0.4λ or I/I sat 16.
We next employ SCEs [Eqs. (5) and (6)] to analyze the coherent transmission through much larger ensembles, for which the full quantum treatment is inaccessible. In Fig. 3 we show how the transmission lineshapes of a 10 × 10 array significantly differs from the Lorentzians of independent atoms. For a single atom the linewidth is power broadened to γ PB (I) = γ √ 1 + I/I sat . In the interacting case, the coherent lineshape is also power broadened but, depending on whether the linewidth of the dominant symmetric collective eigenmode υ at low light intensity [64] is subradiant (υ < γ) or superradiant (υ > γ), it will also be narrower or broader [ Fig. 3(b)], respectively, than γ PB (I). There is no analogous broadening of the incoherent lineshape, however [70]. many-body lineshape also exhibits a dip or "hole burning" on resonance [ Fig. 3(b)]. This dip is analogous to vacuum Rabi splitting [38,39], where the interatomic DD coupling has now taken the role of the cavity coupling, and, while it only occurs for sufficiently high density, it can interestingly still exist even in the fully random ensemble. 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 [3] need not necessarily be restricted to dilute systems. Moreover, 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 [71], with applications to quantum information processing [29], making understanding nonlinear transmission essential. We calculate this for large arrays in Fig. 3(c), and find that the reduction in the extinction as a function of light intensity is considerable -although less prominent with smaller spac- ings. This may ultimately restrict the applications of atomic arrays as highly reflective cooperative mirrors to weak light intensities only.
To conclude, by comparing SCEs and QME, we have identified quantum effects 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 manybody effects reminiscent of single-atom power broadening and vacuum Rabi splitting. The presence of many-body effects even in uniform disks, that are attracting considerable inter- The single atom Hamiltonian H sys, j (in which we have assumed the rotating wave approximation) has the form where ∆ = ω − ω ge is the detuning of the laser frequency ω from the atomic transition frequency ω ge , d ge = d * eg is the dipole matrix element, E + is the positive frequency component of the laser amplitude (given in terms of the electric displacement D + L = 0 E + ), and the raising operator from the ground state |g to the excited state |e ,σ ( j) + = |e j j g|, lowering operatorσ ( j) − = |g j j e|, and excited state population operatorσ ( j) ee = |e j j e| = 1 − |g j j g| are single-atom operators for the jth atom. 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. (S1), are the real and imaginary parts of the dipole radiation kernel G(r): where is the electric field amplitude for an oscillating electric dipole at the origin,n = r/r, and k = 2π/λ for laser wavelength λ. Note that we typically drop the contact interaction term [S1]. Once the full density matrix ρ {r 1 ,...,r N } (t) for a particular set of fixed atomic positions {r 1 , . . . , r N } is known, the one-body ρ ( j) ab ( jth atom), two-body ρ ( j, ) ad;bc ( jth and th atoms), etc., expectation values for this stochastic realization are given by and so forth.
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 ρ ( j) ab (t), ρ ( j, ) ad;bc (t), etc., to spatial correlation functions for the atoms at any given time t: arXiv:1907.07030v2 [quant-ph] 17 Jul 2019 and so forth for higher-order correlations, where the field operatorsψ † a (r) andψ a (r) create and annihilate atoms in internal state a ∈ {g, e} at position r. The atomic correlation functions for a single realization of fixed atomic positions {r 1 , . . . , r N } (as indicated by the subscript) are given in terms of ρ ( j) ab and ρ ( j, ) ad;bc by and it is through solving the coupled dynamics between the light and atoms for each stochastic run and ensembleaveraging over many such realizations that we establish the light-induced spatial correlations between atoms [S1, S2].
In the limit of low light intensity, the overlap between the incident laser field and the eigenvectors v j of the matrix formed by Eq. (S3) [ignoring the contact term in Eq. (S4)] determines the resonant behavior of the atomic ensemble [S3-S5]. Because the matrix is complex symmetric rather than Hermitian, the collective eigenmodes v j are not necessarily orthogonal, resulting, e.g., in asymmetric Fano-like interference resonances, such as between the in-phase and out-of-phase eigenmodes, as shown in Fig. 1(a).
We denote the eigenvalues of Eq. (S3) by ν j + iυ j , where ν j = ω eg − ω j are the shifts of the collective mode resonances from the single-atom resonance frequency and υ j denote the collective radiative half-width at half-maximum (HWHM) linewidths. For υ j > γ (υ j < γ ) the mode is superradiant (subradiant), where γ is the independent-atom linewidth.
For a single, isolated atom, beyond the low light intensity, the solution to Eq. (S1) (i.e., the optical Bloch equations) in the steady state is where the intensity I = 2 0 c|E + | 2 , the saturation intensity is given by I sat = 4π 2 γc/3λ 3 , and the linewidth γ of the terms in Eq. (S12) and Eq. (S13) experiences a power broadening γ PB = γ √ 1 + I/I sat .
We compare the full quantum solution of QME [Eq. (S1)] with the semiclassical equations (SCEs) for the single-body terms ρ ( j) ab based on the factorization ρ ( j, ) ad;bc (4) and (5) in the main text], which neglects quantum fluctuations. In terms of the stochastic sampling procedure, we express this semiclassical factorization as where the [1 − δ(r − r )] term is necessary to exclude the case where the annihilation operators refer twice to the same atom. Despite the factorization of the internal atomic correlation functions, we generally have ψ † 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 (2D) array with precisely one atom per site; and (ii) a random, uni-form 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 [S3], 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 (S15) where the standard deviations σ x , σ y , σ z quantify the spatial confinement of the trapped atoms in all three directions.

SII. 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 whereP is the atomic polarization. To analyze the different contributions in the scattered light, we write it asÊ + d = Ê + d + δÊ + d , where δÊ + d denotes the fluctuations. We then obtain (S17) 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. (S17) 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. (S17) 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. (S16), the coherently scattered light consists of expec-tation values and there is no difference between the quantum and semiclassical coherent scattering. Hence, any difference between the quantum and semiclassical coherent scattering for a manyatom ensemble is due solely to many-body quantum effects. The incoherent contribution in Eq. (S18) is more subtle, as means the incoherently scattered light from a single atom yields (S21) where we have usedσ +σ− =σ ee . In the semiclassical approximation, where the quantum fluctuations are ignored, one then replacesσ + by σ + in Eq. (S20) [S6], such that and the incoherently scattered light intensity in Eq. (S21) 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. (S20) now becomes where, as in Eq. (S16), [G(r−R)] * acts onP − (R) and likewise G(r − R ) onP + (R ). When calculating the full quantum solution, the correlation functions are evaluated using the solution to QME [Eq. (S1)] and by ensemble-averaging over many realizations of atomic positions. However, we can also introduce the many-body version of the single-atom semiclassical approximation [Eq. (S22)] to light scattering: substituting this back into Eq. (S23) to give the semiclassical scattered field Deriving the semiclassical scattered light in Eq. (S25) corresponds to a systematic way of neglecting all quantum fluc-tuations when the atomic response is first calculated from SCEs [Eqs. (4) and (5)]. Hence, comparing the scattered light of Eq. (S25) with the equivalent full quantum solution of Eq. (S23) 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. (S21)] 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. (S23) in the normal order. This yields for the expectation term on the right hand side of Eq. (S23) (for both fermionic and bosonic atoms) Substituting this into Eq. (S23) and using the semiclassical factorization approximation of Eq. (S14) we obtain where denotes a double integral over all {R, R } excluding R = R . The difference between this augmented (semiclassical plus single-atom quantum) expression and the semiclassical expression of Eq. (S25) in the scattered intensity is effectively the contributions of the single atom incoherent (quantum) scattering from Eq. (S21) integrated over the extent of the sample: 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. In this work we consider coherently and incoherently transmitted light and calculate them through a disk of crosssectional 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 fibre, although, for simplicity, we ignore any explicit refocussing or fibre coupling. The transmitted light therefore has the form Note that because of the double integral over S and S the expectation term is now a function of r and r , although substituting r 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 ≡ − log T coh ), we substitute the first four terms on the right hand side of Eq. (S17) into Eq. (S29). 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 Eq. (S29) with Eq. (S18). Evaluating Eq. (S18) using Eqs. (S23) and (S26), along with the solutions to QME, results in the quantum incoherent transmission T QM inc . Using instead the solutions to SCEs and either Eq. (S25) or Eq. (S27), 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.

SIII. FLUCTUATIONS IN INCOHERENT LIGHT
In Figs. 1 and 2 of the main section we analyzed the incoherent transmission and found it to be well approximated by T SAQ inc , with T SC inc providing a negligible contribution at high intensities. In

SIV. SCATTERING FROM A SINGLE ATOM
To illustrate how many-atom collective response affects transmitted light, we show in Fig. S2 the coherent and incoher-ent transmission for a 2 × 2 array, incoherent transmission for a 10 × 10 array, and for a single atom. For a single atom, the coherent scattering (a-c) is purely classical since T QM coh = T SC coh [Eq. (S19)], while the incoherent scattering (df) is purely quantum since from Eqs. (S21) and (S22) we get T QM inc (T SC inc = 0). For the 2 × 2 array we plot T QM coh (a-c) and T QM inc (d-f); compare these with T SC coh and T SAQ inc of Fig. 1 that use the same parameters as Fig. S2(a-f) For the coherent transmission [ Fig. S2(a-c)], the single atom lineshape is given by a single resonance with powerbroadened linewidth γ PB = γ √ 1 + I/I sat [Eq. (S12)]. The many-atom lineshapes, however, exhibit clear qualitative differences, including multiple resonances and modified powerbroadened linewidths, clearly indicating the effects of the sample geometry and light-mediated interactions.
The incoherent transmission through a 10 × 10 array [ Fig. S2(g-i)] is calculated using the semiclassical model incorporating the single-atom quantum description, T SAQ inc , using the same parameters as the coherent transmission in Fig. 3(a,b) of the main text. Similarly to the full quantum 2 × 2 case in (d-f), the lineshape approaches the single atom quantum lineshape with HWHM 5.5γ.

SV.A. Paraxial Gaussian beam
The amplitude of a paraxial Gaussian laser beam in the absence of atoms propagating in the z direction and focused at z = 0 has the form E + (r) = E + w 0 w e ikz e ikρ 2 z /2R c e −iζ(z) e −ρ 2 z /w 2ˆ , where E + is the maximum amplitude, w = w 0 (1 + z 2 /z 2 R ) 1/2 the beam radius, w 0 the beam waist, z R = πw 2 0 /λ the Rayleigh range, ρ z = (x 2 + y 2 ) 1/2 , R c = z + z 2 R /z the beam curvature, ζ(z) = arctan(z/z R ) the Gouy phase, andˆ the unit polarization vector. In every figure except for Fig. 3(c), the beam waist w 0 = 10λ is sufficiently large for the paraxial model to be a good approximation to the true beam propagation.

SV.B. Vector Gaussian beam
In Fig. 3(c) of the main text we consider beam waists of w 0 ≤ 3λ, at which point the vector nature of the light must be correctly accounted for. This is carried out numerically with the method used in Refs. [S7-S9]. We consider a field E + e −ρ 2 L /w 2 Lˆ + with a Gaussian profile incident on a lens at position z = − f , where f is the focal length of the lens, ρ L and w L are respectively the radial position and beam radius at the lens, andˆ ± = ∓(x ± iŷ)/ √ 2 is a circular polarization unit vector. Immediately after passing through the (ideal) lens, the field has the form 1 (a-c) The optical depth of the coherent transmission calculated with QME, OD QM coh = − log T QM coh and the incoherent transmission calculated using (d-f) QME T QM inc and (g-i) the semiclassical model with single-atom quantum description T SAQ inc . We consider a 2 × 2 (a-f) and 10 × 10 (g-i) array with fixed positions, using parameters, respectively, from Figs. 1 and 3(a,b) where φ = tan −1 (y/x), and θ = tan −1 (ρ L / f ) is the angle between the −z axis and a point on the lens. To calculate the field propagation, it is helpful to decompose the field into an orthogonal set of modes: E + = µ κ µ E + µ , where µ = (k t , s, m), k t = k 2 − k 2 z is the transverse wavevector component, s = ±1 is the helicity, m is an angular momentum index, and the expansion coefficients are given by with the J m describing mth order Bessel functions. The field, taken at a distance z from the lens focus (located at the origin), is then given in terms of the decomposed σ ± and z polarization components by E + + (ρ z , φ, z) =E + s=±1 k 0 dk t sk + k z 4πk J 0 (k t ρ z ) e ik z (z+ f ) κ µ , E + − (ρ z , φ, z) =E + s=±1 k 0 dk t sk − k z 4πk J 2 (k t ρ z ) e ik z (z+ f ) e 2iφ κ µ ,