Propagation of light in cold emitter ensembles with quantum position correlations due to static long-range dipolar interactions

We analyze the scattering of light from dipolar emitters whose disordered positions exhibit correlations induced by static, long-range dipole-dipole interactions. The quantum-mechanical position correlations are calculated for zero temperature bosonic atoms or molecules using variational and diffusion quantum Monte Carlo methods. For stationary atoms in dense ensembles in the limit of low light intensity, the simulations yield solutions for the optical responses to all orders of position correlation functions that involve electronic ground and excited states. We calculate how coherent and incoherent scattering, collective linewidths, line shifts, and eigenmodes, and disorder-induced excitation localization are influenced by the static interactions and the density. We find that dominantly repulsive static interactions in strongly confined oblate and prolate traps introduce short-range ordering among the dipoles which curtails large fluctuations in the light-mediated resonant dipole-dipole interactions. This typically results in an increase in coherent reflection and optical depth, accompanied by reduced incoherent scattering. The presence of static dipolar interactions permits the highly selective excitation of subradiant eigenmodes in dense clouds. This effect becomes even more pronounced in a prolate trap, where the resonances narrow below the natural linewidth. When the static dipolar interactions affect the optical transition frequencies, the ensemble exhibits inhomogeneous broadening due to the nonuniformly experienced static dipolar interactions that suppress cooperative effects, but we argue that, e.g., for Dy atoms such inhomogeneous broadening is negligible.

Disordered media of resonant light emitters provide rich mesoscopic systems [24][25][26] where light can mediate strong interactions.Cold and dense atomic ensembles in the limit of low light intensity (LLI) represent systems of dipolar emitters where light can induce strong position-dependent correlations, even though each individual atom's response to a coherent light field behaves like a linear classical oscillator [27][28][29].The light-mediated resonant DD interactions are highly sensitive to atomic positions, leading to radiative excitations in each atom that are correlated to the positions of every other atom in the sample.This granularity of atoms, which determines their optical response, cannot be captured by standard electrodynamics of continuous media, where light-induced interactions between atoms are treated in an averaged sense [29][30][31][32].In the limit where the atoms are considered stationary during imaging, the collective optical response of N atoms is dependent upon the position correlation functions ρ j (r 1 , . . ., r j ), for j = 1, . . ., N, that are determined before the light interacts with the sample [28].When atomic positions are independent and uncorrelated, all these initial correlations factorize ρ j (r 1 , . . ., r j ) = ρ 1 (r 1 ) . . .ρ 1 (r j ), and solving for the optical response involves calculating the radiative excitations of the atoms that can, nevertheless, still be correlated in a disordered system [33,34].However, highly non-trivial position correlation functions between electronic ground-state atoms ρ j (r 1 , . . ., r j ) can exist due to quantum statistics and atomic interactions.Aside from the fundamental interest in cooperative responses in such systems, the scattered light then also directly conveys information about these correlations, serving as a diagnostic method.
Here we solve the collective optical responses of stationary dipolar emitters experiencing static long-range DD interparticle interactions and discuss experimental consequences of the findings.The emitters can represent atoms or molecules, provided that radiative coupling between molecular vibrational states can be neglected [35].The positions of dipoles, determined in the absence of light at zero temperature, are sampled using quantum Monte Carlo methods [36,37] which approximately generate the position correlation functions ρ j (r 1 , . . ., r j ), encompassing all orders j = 1, . . ., N. We consider strongly confined traps that take the shape of oblate (pancake-shaped) and prolate (cigar-shaped) geometries, where the static DD interaction is dominantly repulsive.We focus on the interaction regimes that are sufficiently weak, such that the density distributions are not crystallized.The configurations of dipoles where the positions do not fluctuate independently dramatically influence the optical response in situations where the system is homogeneously broadened, as previously also demonstrated in simulations of light propagation using the Metropolis algorithm in the presence of nontrivial correlations ρ j (r 1 , . . ., r j ) for quantum degenerate ideal fermionic atoms in a one-dimensional optical waveguide [38].The effect of static DD interactions on light scattering has recently been studied in Ref. [39] without including position correlations ρ j (r 1 , . . ., r j ) of the atoms in the analysis.
The strength of recurrent scattering in cooperative optical responses depends on the interparticle separation in units of arXiv:2310.16158v3[cond-mat.quant-gas]11 Jan 2024 the resonance wavenumber of light.Increased density leads to more pronounced deviations from the single-atom Lorentzian lineshape.To characterize the impact of the static DD interactions, we systematically vary the strength of these interactions while approximately maintaining a constant peak density.Our observations reveal that, in both prolate and oblate traps, repulsive static interactions lead to short-range ordering among the dipoles, which, in turn, curtails fluctuations in the light-induced resonant DD interactions, particularly as these interactions become very large at short interparticle separations.This phenomenon is identified in increased coherent reflection and optical depth that are accompanied by reduced incoherent scattering.In an oblate trap, coherent transmission and reflection resonances narrow at low density but broaden at high density.In a prolate trap, the scattered light resonance can narrow even below the natural linewidth.The collective resonance shifts are substantially smaller than those predicted by the collective Lamb shift in continuous medium electrodynamics, underscoring the violation of the standard optics in the system [29,30,40,41].Intriguingly, especially in a prolate trap, the presence of the static DD interactions enables better targeted excitation of subradiant eigenmodes at high densities, although the linewidth of the eigenmodes can exhibit significant variation between different realizations due to the position fluctuations in the ensemble.Additionally, we find that the main difference between the two-level and the isotropic J = 0 → J ′ = 1 transitions is the emergence of resonances where the orientations of the dipoles in an oblate trap, in the latter case, are parallel to the normal of the trap plane.
In disordered dense samples, the excitations can become very localized leading to high concentrations of energy.We analyze the excitations peaks and their widths in individual realizations and find that the static interactions enhance the localization peak strengths, providing control and manipulation of optical fields on a subwavelength scale.We also study the effect of the static DD interactions on optical transition frequencies.When these interactions start substantially influencing the resonances, atoms become inhomogeneously broadened due to the effect of the atoms experiencing the static DD interactions nonuniformly.This broadens the resonances and reduces cooperativity.However, the simulations indicate that such effects for Dy atoms are negligible.
The article is organized as follows: Section II provides the theoretical background for describing static interactions and light-matter coupling, while Sec.III gives an overview of stochastic electrodynamics and the quantum Monte Carlo methods.The optical responses in an oblate trap, when the strength of the static interactions is varied for a constant peak density, are in Sec.IV A, the differences between twolevel and isotropic transitions in Sec.IV B, inhomogeneous broadening in Sec.IV C, and the results in a prolate trap in Sec.IV D. The localization of excitations is presented in Sec.IV E and the optical responses for variable densities in Sec.IV F before concluding remarks in Sec.V.

A. Static long-range dipole-dipole interactions
We assume particles, henceforth referred to as atoms, with the mass M that are harmonically trapped and experience long-range DD interactions where the trapping potential with the frequency ω j , is associated with the characteristic length scale ℓ j = [ℏ/(Mω j )] 1/2 .We consider the atoms confined in an oblate or prolate trap (Fig. 1) with all the static (e.g.magnetic) dipoles µ j oriented in the same direction.The static DD interaction potential V dd between the atoms, which is independent of the light-induced radiative optical DD coupling between the atoms, is then given by where r ℓ j = r j − r ℓ is the vector joining the atoms j and ℓ and rℓ j = r ℓ j /|r ℓ j |.For magnetic dipoles, we have C = µ 0 .The DD interaction can be characterized by the interaction length The DD interaction diverges at small interatomic separations and, following Ref.[12], we introduce an additional repulsive short-range interaction potential, similar to the Lennard-Jones potential, We maintain in the simulations fixed ratios c 6 = 0.0271R 6 dip E dip and c 12 = 4.47 ] where V dd is dominantly repulsive.

B. Collective atom-light interactions
We solve the optical response of the atoms in the limit of LLI where the individual atoms respond to light linearly and behave as a set of coupled, driven harmonic oscillators.In this section, we first consider the commonly employed basic formalism [42] for the coupled dynamics between the atoms and light for a fixed set of atomic positions {r 1 , r 2 , . . ., r N }.When the positions fluctuate, these equations are then utilized in stochastic electrodynamics simulations for a disordered system in Sec.III A. The atoms are illuminated by a near monochromatic Gaussian incident light field with the positive frequency component of the amplitude E(r), the wavevector k, and frequency ω = kc = 2πc/λ.The incident field drives a two-level or a |J = 0, m J = 0⟩ → |J ′ = 1, m J ′ = µ⟩ atomic transition and is detuned ∆ ( j) µ = ω − ω ( j) µ from the resonance ω ( j) µ of the level µ in the atom j that may vary between the atoms due to the interactions.The system dynamics is described by the positive frequency components of the atomic polarization amplitudes P of the dipole matrix elements d j = DP ( j) .All observables are expressed in terms of slowly varying field amplitudes and atomic variables by factoring out the rapidly oscillating dominant frequency of the driving light Ee −iωt → E, P ( j) e −iωt → P ( j) [28].
The positive frequency component of the scattered field E (ℓ) s from the atom ℓ at r ℓ is given by where G(r) is the monochromatic dipole radiation kernel [43] describing the radiated field at point r due to an oscillating dipole at the origin with the dipole d where r = |r| and r = r/r.The positive frequency component of the total external field driving the atom j is given by the incident field plus the scattered fields from all the other atoms Thus, the dipole amplitude of each atom depends on the fields from all other atoms in the system, leading to a set of coupled equations for N atoms with a configuration of positions {r 1 , r 2 , . . ., r N }.In the LLI limit of the coherently driven system, these equations determine the polarization of each atom and the population of the excited state can be neglected [33].
The dynamics then corresponds to a damped classical linear oscillator amplitude driven by E ext dP ( j) where G ( jℓ) µν = ê * µ • G(r j − r ℓ )ê ν and the single atom linewidth We solve Eq. ( 9) in the steady state, assuming that the time scale 1/γ is short compared with the time scale of any centerof-mass motion of the atoms.The dynamics of the atomic polarizations can be compactly represented in the matrix form [34] where b 3 j−1+ν = P ( j) ν , f 3 j−1+ν = iDê * ν • E(r j )/ℏ.The offdiagonal elements of the matrix H describe the light-mediated interactions between the atoms and are given by The diagonal elements are iγ, while the diagonal matrix δH contains the detunings ∆ ( j) µ .Finally, the total field can be calculated from the sum of the incident light field and scattered field from all the atoms We describe the response in terms of the collective excitation eigenmodes v j which correspond to the eigenvectors of H [44][45][46], with eigenvalues δ j + iν j .These represent the collective line shift δ j = ω 0 − ω j from the single atom resonance ω 0 (here assumed degenerate) and the collective linewidth ν j .When ν j < γ, excitations are called subradiant, whereas ν j > γ, they are super-radiant [47].We determine the occupation of the eigenmode b from [42,48]

III. STOCHASTIC SIMULATIONS A. Stochastic electrodynamics
We have discussed how the coupled dynamics between atoms at a fixed set of positions {r 1 , r 2 , . . ., r N } and coherent light in the limit of LLI can be solved using Eqs.( 9), (13), and (6).When atomic positions fluctuate, the system becomes a disordered medium where light can establish correlations among the atoms in cold and dense ensembles that significantly modify the optical response, even within the classical regime [29].Here we briefly review how position fluctuations are conveniently described by employing secondquantized atomic field operators for the ground and excited states ψg,e (r) [28], where the label e implicitly incorporates the Zeeman levels of the J = 0 → J ′ = 1 transition.The positive frequency component of the atomic polarization density operator then takes the form P(r) = ψ † g (r) d ge ψe (r), with d ge representing the dipole matrix element between the electronic ground and excited levels.We adopt the convention that the repeated level index e is summed over.The expectation value of P(r) with respect to the fixed atomic positions {r 1 , r 2 , . . ., r N } is related to P ( j) in Eq. ( 9) via [34] ⟨ Solving for the light field from Eqs. ( 13) and ( 6) in the general case of fluctuating positions requires an integral of the radiation kernel and the expectation value ⟨ P(r)⟩.Multiply scattered light establishes correlations between the atoms, leading to a hierarchy of equations of motion for the correlation functions of atomic density and polarization [27,28].Specifically, in the limit of LLI, we introduce normally ordered correlation functions The LLI is indicated by including at most one P factor in each expression and by the fact that the ground-state correlations ρ k are not perturbed by light.Then P k (k = 1, . . ., N) satisfy the dynamics for degenerate, uniform transition frequencies The ground-state correlations ρ p exist before the light enters the sample.P p (r 1 , . . ., r p−1 ; r p ) represents the correlation function for ground state atoms at r 1 , . . ., r p−1 , given that there is polarization at r p .The strong coupling emerges from the third term on the right-hand side, representing repeated exchanges of a photon between the same atoms, giving rise to recurrent scattering processes.The challenge of solving Eq. ( 18) arises from the complex hierarchy of N equations for N atoms.Standard optics of continuous media is recovered by factorizing all correlations [28][29][30][31], i.e., disregarding any light-established correlations between the atoms P k (r 1 , . . ., and considering an initial classical uncorrelated distribution of atoms ρ n = ρ n 1 .At very low densities, ρ/k 3 ≪ 1, the hierarchy can be truncated, e.g., by including only correlations between pairs of atoms, resulting in a closed equation for P 2 .In this case, the optical response also depends on the groundstate pair correlation function ρ 2 (r 1 , r 2 ).The behavior of a resonance linewidth is strongly geometry-dependent but can be integrated for semi-infinite quantum degenerate ensembles, already exhibiting through ρ 2 (r 1 , r 2 ) linewidth broadening for bosonic atoms [27], narrowing for fermionic atoms [49], and the effects of quasiparticle pairing [50,51].
An alternative approach to solving the LLI response of Eq. ( 18) is through a stochastic Langevin-type method.This method considers the dynamics for excitations for a given set of fixed atomic positions {r 1 , r 2 , . . ., r N } while treating the atomic positions as stochastic variables.In each stochastic realization, the positions of the atoms are sampled to match the proper correlation functions ρ p (r 1 , . . ., r p ), for p = 1, . . ., N in the absence of light.Subsequently, the excitations and the scattered light for the specific set of atomic positions {r 1 , r 2 , . . ., r N } are solved using Eqs.( 9) and (13).The expectation values are obtained by ensembleaveraging over many realizations.This stochastic classicalelectrodynamics approach of coherent scattering can be formally shown to converge to an exact solution for stationary atoms at arbitrary densities by reproducing the correct hierarchy of correlation functions for both 1D scalar theory [33] and 3D vector-electrodynamics with a single electronic ground level [34].The probability distribution required to sample the normally-ordered ground-state atom correlation functions ρ p (r 1 , . . ., r p ) is determined by the many-body wave function P(r 1 , r 2 , . . ., r N ) = |Ψ(r 1 , r 2 , . . ., r N )| 2 [33].
For uncorrelated atoms, where ρ n = ρ n 1 , each atom can be sampled independently.Such scenarios include an ideal Bose-Einstein condensate, a Mott-insulator ground state in an optical lattice, and classical atoms.In the case of atoms obeying Fermi-Dirac statistics in a 1D waveguide, the optical response has been solved by stochastically sampling the atomic positions [38].In Ref. [38], the atomic positions are defined by the Slater determinant of the many-body wave function and sampled using the Metropolis algorithm.Additionally, to account for the influence of a hard-core radius of the classical atom distribution, simulations have been conducted in light scattering by excluding a hard-sphere volume around the already existing atoms in the sampling [52].
Here we calculate the optical response of atoms that exhibit quantum-mechanical position correlations due to static longrange DD interactions that are present before the light enters the sample.We consider an atomic gas in a strongly confined oblate or prolate trap at zero temperature where the DD interactions are dominantly repulsive.These repulsive static DD interactions inhibit small interatomic separations, therefore suppressing light-mediated resonant DD interactions through increased interatomic spacing.While this effect draws some parallels to Fermi-Dirac statistics [38,49], the DD repulsion between the atoms quickly becomes more substantial and is more long-range, which obscures quantum-statistical characteristics of the atoms in the optical response.Although the simulations are performed with bosonic atoms, any signatures of Bose-Einstein statistics in the scattering are rapidly washed out by the DD interactions.The ensemble can represent a strongly correlated quantum many-body system with long-range static DD interactions where the position of every atom affects the positions of every other atom.The correlated atomic positions in the calculations are therefore sampled using quantum Monte Carlo methods, as described in the following section.

B. Quantum Monte Carlo sampling
Solving the stochastic electrodynamics simulations of the optical response, described in Sec.III A, requires sampling the atomic positions in a dipolar gas that are distributed according to the square modulus of the ground-state wave function at zero temperature.We have used the variational and diffusion quantum Monte Carlo (VMC and DMC) methods to find approximate ground-state wave functions and stochastic realizations of atomic positions to solve for the radiative excitations in Eq. ( 9) and scattered light in Eqs. ( 13) and ( 6) for each given set of fixed positions.
In the VMC method, quantum mechanical expectation values are evaluated using a trial many-body wave function that is an explicit function of interparticle distances.The Metropolis algorithm is used to sample position realizations from the square modulus of the wave function for Eq. ( 9), estimators of the energy and other expectation values.Free parameters in the trial wave function are obtained by minimizing the energy expectation value [53].
In the DMC method, drift, diffusion and branching/dying processes governed by the many-body Schrödinger equation in imaginary time are simulated in order to project out the ground-state component of a trial wave function [36].The product of the trial wave function and the solution of the imaginary-time Schrödinger equation is represented as the ensemble average of a discrete population of walkers (weighted delta functions), and the Green's function for the imaginary-time Schrödinger equation is treated as a transition-probability density for walkers over discrete time steps.After equilibration, the walkers are distributed as the product of the trial wave function and its ground-state components, and walkers provide atomic configurations for Eq. ( 9) and ensemble-averaging expectation values.For a bosonic gas the ground-state wave function is nodeless and hence there are in principle no uncontrolled approximations in the DMC method.In practice, the trial wave function must have a sufficiently large overlap with the ground-state wave function that the algorithm can be equilibrated on a tractable timescale.
The DMC algorithm generates atomic configurations distributed as the product of the trial wave function and its ground-state component.The error in the distribution of atomic configurations is therefore first order in the error in the trial wave function.Our trial wave functions were of form where the orbitals ϕ(r i ) were Gaussian functions that could be (initially) centered in the middle of the trap obtained by a brute-force minimization of the potential energy.The position and the width of each Gaussian orbital in each Cartesian direction were treated as optimizable parameters.The Jastrow exponent J contained polynomial two-and three-body terms [54].In addition, we used two-body Jastrow terms designed to impose physically appropriate behavior on the wave function at short range.If the interaction between atoms were of the isotropic form d 2 /r 3 , where d is a constant, then the Jastrow exponent would have to contain a pairwise term u d (r) = − 8d 2 µ/r to ensure that the local-energy contribution from a coalescing pair of atoms, E L2 = −∇ 2 Ψ/(2µΨ)+d 2 /r 3 , diverges more slowly than 1/r 3 at short range, where µ is the reduced mass of a pair of atoms.This negative, divergent two-body Jastrow term u d (r) has the effect of making the wave function go rapidly to zero at coalescence points, and we have continued to use this term even in the presence of anisotropic dipolar interactions.Similarly, two-body Jastrow terms going as −1/r 5 were used to impose the exact short-range behavior on the wave function in the presence of a repulsive r −12 interaction in our calculations using Lennard-Jones potentials.All our VMC and DMC calculations were performed using the casino software [37].

C. Scattered light
Since we consider here the LLI limit where each atom responds to light classically, the scattered light from the atomic ensemble is directly evaluated from the stochastic simulation data by ensemble-averaging over many realizations [55].We then have for the transmission T and reflection R where E and E * denote the positive and negative frequency components, respectively, and the intensity I = 2cϵ 0 ⟨E * • E⟩.
The lens regions centered on the optical axis in the forward and backward directions are A and A ′ , respectively.The scattered field consists of a mean field ⟨E s ⟩ and fluctuations δE s = E s − ⟨E s ⟩ [56].The coherent scattering originates from the mean field part, while the incoherent contribution is due to the position fluctuations of the atoms.Then where the last term includes both the coherent (the first term) and incoherent (the second term) contributions.The coherent scattering is mostly directed in a narrow cone along the incident field direction, while the incoherently scattered light is less directed, and we vary the numerical aperture (NA) of the lenses in these two cases.The coherent transmission is frequently expressed as the optical depth OD coh = − ln T coh .

IV. OPTICAL RESPONSE
We first consider the optical response of weakly excited, independently sampled trapped atoms.This corresponds to a cold disordered atomic ensemble where the atomic positions are not correlated, but light can establish correlations between the excitations of different atoms at small interatomic separations due to random position fluctuations.Variations of related scenarios have been investigated theoretically and several experiments have reached at least close to the regime where a random ensemble of atoms could show such correlations [29-32, 40, 41, 46, 57-78].The optical response in an oblate trap is illustrated in Fig. 2, while the atom number is varied for the J = 0 → J ′ = 1 transition.Even though the atoms are noninteracting before the light enters the sample, at increasing densities due to the light-mediated interactions the resonances are broadened and the lineshapes increasingly deviate from the independent-atom Lorentzian profile, exhibiting multiple peak resonances and asymmetric profiles.Coherent transmis-sion and reflection also show a pronounced density-dependent resonance shift.The dependence of the correlated response on the peak 2D density per wavenumber squared ρ2D /k 2 is clearly visible in Fig. 2. For ρ2D /k 2 ≪ 1, the response is close to the independent-atom Lorentzian, while cooperative recurrent scattering dominates at ρ2D /k 2 ∼ 1.
The effect of the static DD interactions on the atomic density profile is immediately obvious in Fig. 3 where we show the numerically calculated (using Monte Carlo simulations, Sec.III B) ground-state density profile of N = 100 atoms at zero temperature in a symmetric oblate trap.Due to the tight confinement of the atoms along the orientation of the dipoles in the z direction, the lines joining the atoms are almost perpendicular to the dipoles and the interactions are repulsive, as shown in the interatomic dipolar potential in Fig. 3(a).This results in a well-known effect of increasing cloud radii and flattening density profiles.In the simulations of the optical response, we accommodate these effects by varying the beam focusing (also to adjust the overlap between the incident and coherently scattered light) and by using a sufficiently large lens to collect the light.Aside from the density profiles, the static DD interactions have a much more substantial effect on the optical response due to the change in atomic correlations that alter the light-mediated interactions between the atoms, as we will show in the next sections.
FIG. 3. Dipolar potential between two dipoles perpendicular to the separation and the 2D density profile as a function of the radius for 100 atoms for different interaction lengths R dip [Eq.( 4)] for ℓ x /ℓ z = 100.

A. Varying static dipolar interaction in an oblate trap
The crucial parameter in the cooperative response of independent atoms is the atom separation in the units of 1/k that determines the strength of recurrent scattering and the emergence of light-induced correlations (Sec.III A) [27,30].Increasing densities lead to deviations from the continuous media electrodynamics due to these correlations.To isolate the effect of static DD interactions and the resulting position distributions on the optical response, we therefore consider atomic ensembles with constant ρ2D /k 2 while varying the dipolar length R dip [Eq.( 4)].We express the interaction strength as the ratio between R dip and the average separation between atoms in the xy plane at the trap center R dip ρ1/2 2D (in a prolate trap we use R dip ρ1D ).We take N = 200 atoms and adjust the peak 2D density ρ2D by changing the trapping frequencies.The transmission and reflection are shown in Fig. 4 for two-level atoms for linear polarization in the trap plane at the densities ρ2D /k 2 ≃ 0.1 and 1, respectively, by varying R dip ρ1/2 2D ≃ 0, 0.017, 0.15, 1.6.We find an increase in the resonant coherent reflection and optical depth for increasing R dip , as also shown in Fig. 5.This is associated with reduced incoherent scattering (with the exception of the incoherent reflection in the strongest interaction case of ρ2D /k 2 ≃ 0.1 which has particularly high coherent reflection).At high densities, the coherent scattering is already strong for independent atoms R dip = 0.This is not true for the low density case ρ2D /k 2 ≃ 0.1 when the relative increase in the coherent resonant scattering as a function of R dip is more noticeable at weak interaction strengths.At ρ2D /k 2 ≃ 0.1, the transmission and reflection resonance linewidths narrow with increasing R dip , but at the higher density case ( ρ2D /k 2 ≃ 1) there is broadening (Fig. 5).At ρ2D /k 2 ≃ 1, the resonance is shifted as a function of R dip , whereas the interaction-dependent shift at the low density case is absent.At ρ2D /k 2 ≃ 1, the lineshape is substantially deviated from the independent-atom Lorentzian and this effect is magnified by the static interactions.The incoherent scattering exhibits a very broad resonance.
The effect of the static DD interactions on the lightmediated interactions between the atoms at a given density is due to position correlations.Increasing the repulsive interactions modify the distribution of interatomic separations introducing short-range ordering of the atoms, as shown in the probability distributions of the nearest-neighbor and all-atom separations in Figs. 6 and 7.As R dip and repulsion increase, the atoms are prevented from coming close to each other.This suppression of small interatomic separations prevents the light-mediated DD interactions, which scale as ∝ 1/r 3 at small separations, becoming very large.In disordered ensembles this reduces the fluctuations of light-induced resonance shifts between any atom pairs, resulting in reduced inhomogeneous broadening of resonance frequencies due to the resonance shifts and reduced incoherent scattering.In Fig. 7, the most pronounced short-range ordering occurs for R dip √ ρ2D ≃ 1.6 when the atoms are unable to approach each other, representing the most dramatic lineshape deformation in Fig. 4. We additionally find in Figs. 6 and 7 that the large separations are suppressed due to increased trap frequencies that maintain the constant peak density value.

Eigenmodes
The optical response of the atoms can be analyzed in terms of the collective excitation eigenmodes.In Fig. 8, we show the normalized histogram distribution of the eigenmodes as a function of their collective resonance linewidths ν and line shifts −δ for ρ2D /k 2 ≃ 0.1 and 1 (the minus sign is used to align the plots with the laser detuning in lineshapes).This corresponds to the optical responses shown in Fig. 4 for the independent-atom case R dip = 0 and for the most strongly interacting case R dip √ ρ2D ≃ 1.6.Owing to the relatively low density of the ρ2D /k 2 ≃ 0.1 case, the eigenmodes for R dip = 0 are strongly peaked around the single atom resonance and linewidth.The static DD interactions cause the highly occupied region to spread, extending towards blue-detuned superradiance and red-detuned subradiance.This peak region is also shifted consistently with the collective resonance shift in Fig. 4.
With increasing density, the size of the central region generally increases both in width in resonance and in linewidth.The distribution also forms long arms that extend into regions of super-radiant and subradiant modes.At ρ2D /k 2 ≃ 1 in Fig. 8, the mode density is high far from the single-atom resonance even for independent atoms.The distribution becomes highly asymmetric at high density with pronounced blue-detuned subradiant eigenmodes.The static interactions further magnify these modes and also produce very super-radiant reddetuned modes.These changes in eigenmode distributions correspond to resonance broadening of the lineshapes with increasing density in Fig. 4.
We also calculate which eigenmodes are occupied at specific laser frequencies in Fig. 4.Here Eq. ( 14) is used to produce a histogram plot of the occupied mode probability distribution in Fig. 9. Intriguingly, the presence of the static DD interactions allows for better targeted excitation of subradiant eigenmodes at high atom densities, as shown in terms of increased and more localized occupations in Fig. 9(c), (d) for a subradiant eigenmode resonance at the detuning ∆ ≃ −0.6γ.Although the incident light strongly couples to a single eigenmode already for the case of independently distributed atoms, the coupling is even more selective in the presence of static DD interactions.Similarly, due to the static coupling, it is more difficult to excite the modes off resonance.Although the modes are only excited over a narrow range of frequencies, they still extend over a wide range of linewidths due to the position fluctuations even with strong static interactions.
Selectively exciting subradiant eigenmodes of atoms in oblate trapping geometries, as well as storing photons in such modes, has been actively investigated when the atoms form periodic arrays and individual atoms are trapped at specific spatial locations, see, e.g., Refs.[42,48,[79][80][81][82].Here the driving of subradiant modes is studied instead in disordered ensembles with short-range ordering due to the static DD interactions.In Sec.IV D, we show how these effects are even more pronounced in a prolate trap where small atom numbers also permit better targeting of individual eigenmodes.

B. Isotropic vs two-level transition
The optical response of two-level atoms in Fig. 4 is qualitatively similar to the response of atoms with an isotropic J = 0 → J ′ = 1 transition to positive circular polarization in the trap plane.At low densities, the transmission and reflection lineshapes are very similar (Fig. 10).At the high density case with ρ2D /k 2 ≃ 1, the coherent scattering is still only slightly modified, but more notable deviations in the lineshape appear in the incoherent transmission.In the case of the J = 0 → J ′ = 1 transition, the incoherent scattering has a dominant peak at ∆ ≃ 0 and a smaller peak at ∆ ≃ −2γ, while in the two-level case the ∆ ≃ 0 peak is less pronounced.In the isotropic case, the dipoles can be excited in the direction normal to the trap plane, despite these components not being directly driven by the incident light.This can result in scattering that is not captured by the lenses.To calculate this out-of-plane excitation, we define the expectation values of the magnitudes of the atomic polarization components,  ground-state atom distributions and correlations.Depending on the physical system considered, the static DD coupling can also influence the optical transition frequencies.Since the static DD field experienced by each atom depends on the relative positions of the other atoms, the effect is not uniform throughout the sample.For example, in the case of static magnetic DD interactions where the electronic ground and excited levels of each atom exhibit a magnetic dipole, the atoms can experience level shifts due to ground state -ground state, ground state -excited state, and excited state -excited state interactions which depend on the specific level structure and atom. 1 For simplicity, we consider a general model to demonstrate the effect of inhomogeneous broadening of resonance frequencies that results from the nonuniformly experienced static DD interactions between the atoms.We introduce position dependent transition frequency shifts in individual atoms where the shift in the atom j at r j is caused by the static DD interaction of all the other atoms ℓ at positions r ℓ and D denotes the effect of the static field on the transition frequency.Introducing Eq. ( 25) in the simulations broadens the distribution of the atomic transition frequencies in the ensemble and shifts the peak value of the atomic transition frequency as a function of D. In optical responses, the broadening is analogous to inhomogeneous broadening in thermal and other ensembles [29,40,84].The transition frequencies are most strongly shifted for any atom pairs that are very close to each other [see the ground state -ground state coupling in Fig. 3(a)].However, the repulsive force between the atoms reduces the likelihood of very strong shifts; the DD interaction between the electronic ground level atoms inhibits atomic distributions with short interatomic separations.
The nonuniform change in the resonance frequencies broadens the transmission and reflection resonances at different coupling strengths in Fig. 11.To make the different cases comparable, the resonances are shown with respect to the most likely shift in the transition frequency δ in each case.The consequences of these level shifts on Dy atoms are discussed in Sec.V where we argue that their effect on the optical response is negligible.

D. Varying static dipolar interaction in a prolate trap
We consider N = 10 atoms trapped in a prolate trap at peak densities ρ1D /k ≃ 0.1 and 1.The static dipoles again are all parallel and repulsive, oriented perpendicular to the long trap axis.The atoms, with the J = 0 → J ′ = 1 transition, are illuminated by a Gaussian beam, with positive circular polarization, propagating either along the long axis of the trap ('on-axis') or perpendicular to it ('off-axis').The on-axis scattered light power and the atomic pair distributions are shown in Fig. 12.At low atom densities, the lineshape is close to the Lorentzian with the resonance near the single atom resonance, while at high densities it is asymmetric.The coherent scattering resonance is shifted towards negative detuning.Similarly to the oblate case in Sec.IV A, we find the resonance narrowing with increasing R dip (Fig. 13), increasing coherent scattering, and decreasing incoherent scattering, but now the HWHM resonances are very narrow and below the linewidth of an isolated atom.These changes due to R dip originate from the short-range ordering of the atoms that suppresses the fluctuations of the light-mediated DD interactions.
Due to small atom numbers and eigenmodes, it is easier to identify individual modes at high densities in a prolate trap than in an oblate one.The possibility of highly selected targeting of the eigenmodes as a result of static DD interactions, highlighted in Sec.IV A, becomes even more pronounced in a prolate trap, as shown in Fig. 14(e), (f) for ∆ = −0.9γ.We also find that the most occupied eigenmode in most stochastic realizations for the system of Fig. 14(f) has close to 0.9 normalized occupation.The average occupation of the most occupied eigenmode over many realizations is about 0.66 (with the standard deviation of 0.16), as a result of some realizations having two eigenmodes more highly occupied.The linewidth of the excited subradiant mode for N = 10 atoms in this case considerably varies between stochastic realizations, but some subradiant modes with a broader resonance in Fig. 14(f) have a smaller wavenumber than those with a narrower resonance.
The scattering in the off-axis case is markedly different to the on-axis case (Fig. 15).While the coherent scattering is again strengthened by the DD interactions, the lineshape for the off-axis scattering becomes notably deformed, with double and triple peaks appearing in the coherent and incoherent scattering, respectively.The extra incoherent peak indicates the effect of the multilevel structure of the atoms, as the reso- nance is not captured by the lens for coherent scattering.
At the single-atom resonance, the static DD interactions induce subradiant excitations in Fig. 16, but the dominant excitations still exhibit the single-atom linewidth as in the independent-atom case.This unchanged dominant occupation across the interaction strengths explains the unshifted central peak in Fig. 15.The potential for targeted excitation of subradiant eigenmodes due to the static DD interactions is again shown in Fig. 16(c), (d) at large detunings.However, the subradiant modes in the both cases co-exist with superradiant eigenmodes.These super-radiant modes are visible in Fig. 16(c), (d) as distributions for ν > γ that extend over a wide range of resonance frequencies and can be excited even off resonance.It is these super-radiant eigenmodes that are also responsible for the additional resonances in the incoherent lineshape in Fig. 15.Owing to their broad resonance linewidths, these modes show up in the lineshape profile.

E. Localization of atomic polarization
One of the immediate consequences of cooperative emitter responses is the resulting intricate interplay between the collective excitation eigenmodes and disorder in emitter positions.This can dramatically affect the near-field landscape of the optical response, resulting in the localization of eigenmodes and highly concentrated excitation energies [85,86].In optics, this can be exploited, e.g., in achieving strong coupling between the excitation modes and a quantum emitter, thereby modifying its decay rate [87].A strongly localized excitation can effectively drive quantum emitters by acting as a cavity with its quality factor determined by the collective linewidth of the mode.We find that the near-field excitation landscape of the atoms exhibits strong localization of excitations on a subwavelength scale that depends on the atom density and the static interaction strength.The localization rapidly increases with the atom density, independently of the static interaction.Examples of single stochastic runs of atomic polarization density profiles are shown in Fig. 17.In the low density case, the profile is closer to the Gaussian, while the localized excitations are visible at high densities.As the static interaction is increased, the value of the peak excitation increases (Fig. 18).This is more pronounced for the most localized 2%.

F. Varying the atom density in an oblate trap
Instead of fixing the peak atom density, we now study the optical response for fixed static DD interaction length in terms of the resonance wavenumber of light R dip k [Eq.( 4)].This allows to investigate how light transmission and reflection vary with the density by changing the trap frequency ℓ x k for constant ℓ x /ℓ z = 25 (Fig. 19).The setup is similar to the one studied in Sec.IV D (see Fig. 1) and we consider the independent atoms R dip = 0 and static dipoles with R dip k ≃ 0.1.The independent-atom case behaves qualitatively similarly to the previous example where we only varied the atom number.The lineshapes at the peak densities ρ2D /k 2 ≃ 1.4 and 3.3 in Fig. 19 are asymmetric in all the cases, displaying an increased optical depth for blue detuning and enhanced incoherent scattering for red detuning.The static interactions sig- nificantly increase the red-detuned incoherent scattering.The density-dependent resonance broadening in Fig. 20 becomes more dramatic with the static DD interactions.The peak value of coherent reflection first significantly increases with density and then at high densities start decreasing again.This may be due to the eigenmode resonances becoming spectrally more distinguishable, as illustrated for the eigenmode occupations in Fig. 20.At high density ρ2D /k 2 ≃ 3.3, the occupation is prominent only around modes for which the laser frequency is resonant.There is little occupation of modes away from this resonance.This selectivity is responsible for the broad transmission resonances at high densities, as different eigenmodes can be excited at different frequencies and the mode resonances extend over a wide range of frequencies.

V. CONCLUDING REMARKS
We successfully solved a challenging hierarchy of N equations (18) for the correlation functions of N atoms.These equations encompass correlations among atoms in their electronic ground states as well as those involving both ground and excited states.This was made possible using stochastic electrodynamics simulations of coupled radiative dipole excitations where the positions of the dipoles are correlated by static repulsive DD interactions and sampled using quan- tum Monte Carlo methods.For the bosonic fluid-like states considered in this paper, ergodicity is not expected to pose a challenge when seeking the ground state using diffusion quantum Monte Carlo simulation, with the accuracy of the sampled distributions scaling linearly with the error in the trial wave function.As long as the sampling of positions is precise, the stochastic electrodynamics simulations of coherently scattered light converge to an exact solution for stationary, laserdriven atoms at arbitrary densities in the limit of LLI [33,34].The methodology can be extended to include also stronger static interactions and other strongly correlated ensembles.
Our main findings showed how the repulsive static interactions lead to short-range ordering among the dipoles, which, in turn, curtails fluctuations in the light-induced resonant DD interactions (for an oblate trap in Sec.IV A and in a prolate trap in Sec.IV D).This phenomenon affects the resonance widths and shifts (Secs.IV A and IV D).Furthermore, it is identified in increased coherent reflection and optical depth that are accompanied by reduced incoherent scattering.The effects of the static interactions on the optical response can be analyzed in terms of the collective excitation eigenmodes (Secs.IV A 1 and IV D) and we found, e.g., that the presence of the static DD interactions enables much better targeted excitation of subradiant eigenmodes at high densities especially in a prolate trap, despite disordered atom distributions.For an isotropic level structure, the excited dipoles can exhibit nonnegligible orientation even along the normal of an oblate trap (Sec.IV B).Intriguingly, the static interactions can also enhance the peak strengths of the disorder-induced excitation energy localization, providing control and manipulation of optical fields on a subwavelength scale (Sec.IV E).
Optical responses in the presence of static DD interactions could be investigated, e.g., using atoms or polar molecules.The repulsive static DD interactions in a prolate and oblate trap can suppress inelastic losses in both systems even at high densities [1].Dy atoms have a large magnetic moment µ ≃ 10µ B , where µ B is the Bohr magneton [4].For example, the 626nm transition has been used in magneto-optical trapping of 162 Dy [88].For ρ2D /k 2 ≃ 3.3, this gives the average interatomic separation 0.09λ at the center of the trap and R dip √ ρ2D ≃ 0.4, with R dip ≃ 21nm.If the static DD interactions strongly affect the optical transition frequencies, resonant emitters can experience them nonuniformly and become inhomogeneously broadened (Sec.IV C).The strength of the level shifts in Eq. ( 25) generally depends on the level structure, but similar estimates for Dy indicate a much weaker effect than any of the cases shown in Fig. 11, and therefore a negligible influence on the optical response.
For alkali-metal atoms the magnitude of the magnetic dipole |µ| ≲ µ B and R dip is more than two orders of magnitude smaller than for the case of Dy.The simulation results in Fig. 4 then indicate that the static magnetic DD interactions between the atoms could not contribute to observable effects, e.g., in the coherent light scattering experiments of Ref. [71].
Heteronuclear polar molecules possess large electric dipole moments whose strength can be controlled by orienting the molecule by external electric field [16].Also atomic Rydberg excitations under conditions of electromagnetically-induced transparency can be used to manipulate collective optical interactions [89,90] and the DD interaction between atoms in Rydberg states can be tuned to a weak-interaction regime [91].Rydberg transitions induce DD interactions resulting in level shifts that determine which atoms engage in resonant scattering of light.Consequently, the positions of the atoms that are resonant with the incoming light can become correlated due to the Rydberg DD interactions.Furthermore, in recent experiments, it has been demonstrated that polar molecules can be utilized in effectively controlling atomic resonances [9].Although the control in Ref. [9] was obtained through Rydberg states, one could envisage a scenario in which an oblate trap of polar molecules exhibiting repulsive static DD interactions is placed on the top of a prolate atom trap.The interactions between these molecules and the atoms could potentially be harnessed to ascertain, through induced atomic level shifts, which atoms can engage in resonant interactions with light.This selective process might ensure that only the atoms located atop the molecules are in resonance.Consequently, the position correlations initially associated with the molecules would be transferred to the optically interacting atoms.

FIG. 1 .
FIG. 1. Schematic illustration of the trapping geometries illuminated by coherent light.(a) An oblate (pancake-shaped) trap in the xy plane (ℓ z ≪ ℓ x = ℓ y ), with the static atomic dipoles aligned perpendicular to the plane along the light propagation direction.An elongated prolate (cigar-shaped) trap (ℓ z ≫ ℓ x = ℓ y ) with the static dipoles aligned perpendicular to the long axis.The incident light in (b) propagates along either the long or short axis of the trap.

FIG. 12 .
FIG. 12. Coherently and incoherently forward-scattered power (in units of incident light intensity/k 2 ) from 10 atoms in a prolate trap (ℓ z /ℓ x = 25) for different static dipole interaction strengths at the peak densities (a), (b) ρ1D /k ≃ 0.1 and (c), (d) ρ1D /k ≃ 1, and (e), (f) the atomic pair distributions.The light propagates parallel to the long trap axis.Normalized (e) nearest-neighbor and (f) all-pair distributions between the atoms.The lens NA 0.8.

FIG. 17 .
FIG. 17.Localization of excitations |P| [in units of DE/(ℏγ)] as a function of the distance from the most excited atom |ϱ − ϱ max | for 200 atoms in an oblate trap at the peak density ρ2D /k 2 ≃ 1 for (a) all data, for (b) 2% of the most localized cases (1000 realizations).The corresponding profile in an individual stochastic realization with R dip √ ρ2D ≃ 0.15 for (c) ρ2D /k 2 ≃ 0.04, (d) 1.