The maximum refractive index of an atomic medium

It is interesting to observe that all optical materials with a positive refractive index have a value of index that is of order unity. Surprisingly, though, a deep understanding of the mechanisms that lead to this universal behavior seems to be lacking. Moreover, this observation is difficult to reconcile with the fact that a single, isolated atom is known to have a giant optical response, as characterized by a resonant scattering cross section that far exceeds its physical size. Here, we theoretically and numerically investigate the evolution of the optical properties of an ensemble of ideal atoms as a function of density, starting from the dilute gas limit, including the effects of multiple scattering and near-field interactions. Interestingly, despite the giant response of an isolated atom, we find that the maximum index does not indefinitely grow with increasing density, but rather reaches a limiting value $n\approx 1.7$. We propose an explanation based upon strong-disorder renormalization group theory, in which the near-field interaction combined with random atomic positions results in an inhomogeneous broadening of atomic resonance frequencies. This mechanism ensures that regardless of the physical atomic density, light at any given frequency only interacts with at most a few near-resonant atoms per cubic wavelength, thus limiting the maximum index attainable. Our work is a promising first step to understand the limits of refractive index from a bottom-up, atomic physics perspective, and also introduces renormalization group as a powerful tool to understand the generally complex problem of multiple scattering of light overall.

One interesting observation is that all the optical materials that we know of, with a positive index of refraction at visible wavelengths, universally have an index of order unity, n ∼ O(1). While we typically utilize materials far from their natural electronic resonances, this even holds true close to resonance [1][2][3][4][5][6][7][8]. Yet, despite the profound implications that an ultra-high index material would have for optical technologies, a deep understanding of the origin of this apparently universal behavior seems to be lacking. Furthermore, this property of real materials is not readily reconciled with the fact that a single, isolated atom exhibits a giant scattering cross-section σ sc ∼ λ 2 0 for photons resonant with an atomic transition of wavelength λ 0 ( Fig. 1-a), which far exceeds both the physical size of the atom or the typical lattice constant of a solid (λ 0 ∼ 1µm for a typical optical transition, compared to the Bohr radius a 0 ∼ 0.1nm).
Specifically, we investigate in detail the optical response of an ideal ensemble of identical, stationary atoms, as a function of density starting from the dilute limit, and well within the regime where the atoms do not interact chemically. In large scale numerics (involving up to ∼ 25000 atoms, about an order of magnitude larger than comparable works [12, 15, 18-21, 34, 36-40]), we find that the maximum index does not indefinitely grow with density, and saturates to a maximum value of n ≈ 1.7, when the typical distance between atoms becomes smaller than the length scale associated with the resonant cross section, i.e. d < λ 0 . Furthermore, we introduce an underlying theory based upon strong-disorder renormalization group (RG), which has been a very successful technique to deal with highly varying interaction strengths in a wide variety of condensed matter systems [43][44][45][46][47][48][49][50]. In the context of our particular problem, the combination of strong near-field (∼ 1/r 3 ) optical interactions and random atomic positions enables one to characterize the optical response of the system in terms of a hierarchy of strongly interacting, nearby atomic pairs. The shifts of the resonance frequencies arising from the near-field interactions then effectively yield an inhomogeneously broadened optical medium, where the amount of broadening linearly scales with density. This implies that light of any given wavelength only interacts with at most ∼ 1 near-resonant atom per reduced cubic wavelength λ 3 0 /(2π) 3 , regardless of the physical atomic density, thus limiting the optical response ( Fig. 1-d). In a dense ensemble with many atoms per cubic wavelength λ 3 0 , the scattering of an incident photon can involve multiple scattering and interference between atoms. c) In conventional theories of macroscopic optical response, the atoms are approximated by a smooth medium, and the index is derived from the product of single-atom polarizability and density. The maximum index n near the atomic resonance then scales with atomic density like n ∼ N λ 3 0 /V . d) In our renormalization group theory, we retain multiple scattering and granularity, showing that the optical properties of the ensemble are determined by a hierarchy of nearby atomic pairs that strongly interact via their near fields. These interactions effectively produce an inhomogeneously broadened ensemble, where the amount of broadening scales with density (with the different colors of atoms representing the different resonance frequencies in the figure). An incident photon of a given frequency thus sees only ∼ 1 near-resonant atom per reduced cubic wavelength to interact with, regardless of atomic density. This results in a maximum index of n ≈ 1.7.

a) b)
Figure 2. Simulated physical system. a) A cylindrical ensemble of randomly distributed atoms (green points) is illuminated by a z-directed Gaussian beam, whose beam w(z) λ0 is represented in orange. The transverse radius of the cylinder is chosen to be much larger than the beam waist, to avoid edge diffraction. b) 3D representation of the forward scattered intensity I(r, ω0) = |E(r, ω0)| 2 (with the value indicated in the colorbar) over a hemispherical surface far from the ensemble (the radius of this hemisphere is 35λ0), given an input resonant Gaussian beam. The intensity is calculated for a single, random atomic configuration. The system parameters used are: beam waist w0 = 3λ0, cylinder radius lcyl. = 7λ0 and thickness d = 2λ0.
Our results are potentially significant on a number of fronts. First, they provide a convincing picture of why typical theories for optical response, based upon a smooth density approximation, fail for dense, near-resonant atomic media, due to the important role of granularity and nearest-neighbor interactions. Furthermore, our results show the promise of a bottom-up approach to understanding the physical limits of refractive index, starting from objects (isolated atoms) whose optical responses are both huge and exquisitely understood. Separately, the existence of a fundamental mechanism that results in inhomogeneous broadening (i.e. dephasing) and saturation of optical properties at high densities, which occurs even for perfect, stationary atoms, should impose fundamental bounds on the maximum densities and minimum sizes of atom-light interfaces needed to realize high-fidelity quantum technologies. Finally, while we focus here on the linear optical response of a dense atomic medium, we believe that the validity of RG is quite general, and can constitute a versatile new tool for the generally challenging problem of multiple scattering in near-resonant disordered media [12, 20, 31, 34-37, 40, 51-53], including in the nonlinear and quantum regimes [54].
This paper is structured as follows. First, we briefly review the theoretical formulation of the multiple scattering problem of atoms or other point-like dipoles, and the standard atomic physics model of refractive index, when atomic granularity and multiple scattering are ignored. We then formulate our large-scale numerical simulations, describing a few implementation details that allow index to be efficiently calculated, and show that the index eventually saturates with increasing density to a maximum value of n ≈ 1.7. We then introduce our RG theory, which highlights the importance of granularity and nearby atomic pairs on the macroscopic optical response, before concluding with an expanded discussion of future interesting directions to investigate.

I. FORMAL THEORY OF MULTIPLE SCATTERING
We consider a minimal system consisting of N identical, stationary two-level atoms. The atoms are assumed to consist of a ground and excited state |g , |e , with frequency difference ω 0 and associated wavelength λ 0 , and which have an optically allowed transition with a dipole matrix element along a fixed axis (sayx), as depicted in Fig. 1a. The excited states of the atoms decay purely radiatively, with a rate of Γ 0 for a single, isolated atom. As we are specifically interested in the linear refractive index, it is sufficient to treat atoms in the limit of classical, polarizable, radiating dipoles. In order to investigate the frequency-dependent index n(ω), we consider that the atoms are driven by a monochromatic, linearly-polarized input beam E in (r, ω) = E in (r, ω)x, whose polarization aligns with the polarizability axis of the atoms. Each atom j acquires a dipole moment d j (ω) = d j (ω)x, as a result of being driven by the total field, which consists of the sum of the incident field and fields re-scattered from other atoms. Formally, the total field can be expressed as [55] (1) Here, the dyadic Green's tensorḠ(r, r j , ω) encodes the field at position r, produced by an oscillating dipole at r j , and in vacuum is given by [55] G(r, r , ω) = k e iρ 4π with ρ ≡ |ρ|≡ k|(r − r )| and k = ω/c. Note thatḠ(r, r , ω) contains both non-radiative, near-field (∼ 1/ρ 3 ) and radiative, far-field (∼ 1/ρ) terms. Then, the induced dipole moment of atom i is given by where the parameter α 0 (ω) defines the polarizability of a single dipole. Although Eq. 1 and Eq. 3 can describe any system of linearly-polarizable point-like dipoles [55], e.g. dielectric nano-particles [56], in our case we focus on the response of non-absorbing, purely radiative atoms, whose resonant cross section σ sc = 3λ 2 0 /(2π) is the maximum set by the unitarity limit [57]. In this context, the atomic polarizability reads α 0 (ω) = −3π/[(∆ + i/2)k 3 0 ], where k 0 = 2π/λ 0 denotes the resonant wavevector, while ∆ ≡ (ω − ω 0 )/Γ 0 represents the dimensionless detuning between the input beam frequency ω and the atomic resonance ω 0 .
While Eq. 1 and Eq. 3 are formally exact, solving a number of equations that explicitly scales with the number of atoms and that depends on the details of atomic positions is not a particularly convenient way to calculate the index or other optical properties. Historically, this fostered the development of simplified theories for the macroscopic response, such as the Drude-Lorentz model [9] or equivalently the Maxwell-Bloch (MB) equations [10], where the discreteness of atoms is replaced by a smooth medium of density N/V ( Fig. 1-c). The resulting index depends on the product of density and single-atom polarizability, where we defined the dimensionless density η ≡ N/(V k 3 0 ). Notably, for an optimum detuning, the maximum real part of the index scales as ∼ √ η.
While the MB equations ignore multiple scattering, the Lorentz-Lorenz (LL) or the equivalent Clausius-Mossotti model is one well-known approach to approximate its effects, still within the smooth density approximation. Given any atom located at r 0 , the model approximates the neighbouring atoms as a smooth dielectric medium with a small spherical exclusion around r 0 [9]. The resulting local field correction gives an index that satisfies the equation (n 2 LL − 1)/(n 2 LL + 2) = (N/V )α 0 (ω)/3 [55]. Plugging in the atomic polarizability, one readily finds that n LL (∆) = n MB (∆ + πη).
Importantly, while the spectrum is shifted, the LL model still produces a maximum index that grows like ∼ √ η.

II. COUPLED-DIPOLE SIMULATIONS
Eq. 1 and Eq. 3 are ubiquitously used to model multiple scattering and interference effects involving a moderate number of point-like scatterers. Here, we briefly introduce some key details of our implementation, which allows us to perform simulations on very high atom number and efficiently extract the index. First, one conceptually straightforward way to extract the complex refractive index of a material would be to take an extended slab of finite thickness d, and investigate the phase shift and attenuation of a quasi-plane-wave incident field upon transmission. We approximately realize such a situation by taking atoms with a fixed density in a cylindrical volume, illuminated by a weakly focused, near-resonant Gaussian beam ( Fig. 2-a). The parameters are chosen such that the beam waist w 0 is small compared to the radius of the system, so that diffraction effects from the edges are negligible. We also avoid very tight focusing w 0 λ 0 , where non-paraxial effects could emerge.
We must also specify a practical definition of index, for a granular system as ours. In particular, since our atoms are purely scattering and have no absorption, it is well-known [20,51,52,58] that for a fixed random spatial configuration, an input as in Fig. 2 produces a complex "speckle" pattern in the outgoing intensity when the system is optically dense, due to multiple scattering and interference, as exemplified in Fig. 2 To isolate the part of the field that possesses a well-defined phase relationship with the incident field from realization to realization, we project Eq. 1 back into the same Gaussian mode as the input, as can be experimentally enforced by recollecting the transmitted light through a single mode fiber. This results in a transmission coefficient t(∆) given by [36,59] where E 0 is the input field amplitude at the beam focus. Here, for convenience, we have defined re-scaled dipole amplitudes c j (∆) = d j (ω)k 3 0 /(3π 0 E 0 ), which satisfy the dimensionless coupled equations In these equations, we define G ij ≡ (3π/k 0 )x ·Ḡ(r i , r j , ω 0 ) ·x and G jj = i/2, which coincides with the single-atom decay rate in units of Γ 0 , while regularizing the divergent self-energy associated with the real part ofḠ. Note that, for simplicity, the Green's functionḠ(r i , r j , ω 0 ) is only evaluated at the atomic resonance frequency, in order to ease the computational cost as the detuning is varied. Ignoring the dispersion ofḠ is an excellent approximation for nearresonant atoms, as the optical dispersion and delay of such a system is dominated by the atomic response itself rather than from the vacuum [60]. Similarly, we approximate the near-resonant input field as E in (r i , ω) E in (r i , ω 0 ). The expression in Eq. 6 represents a useful closed-form definition of the transmission coefficient t(∆), which avoids a numerically expensive point-by-point evaluation of the scattered field E(r, ω), as nominally prescribed by Eq. 1. We can extrapolate the complex index of refraction n(∆) from the relation where the averages are performed over ∼ 10 3 −10 4 sets of random positions, for each fixed density. Unlike in a smooth medium, we have that | t(∆) | 2 = |t(∆)| 2 . Nevertheless, our definition of the index coincides with that often used within atomic physics (e.g. in phase contrast or absorption imaging of a Bose-Einstein condensate [61,62]). In the Appendix B, we discuss the independence of the calculated index from the thickness d, which is implicitly assumed in Eq. 8. Alternatively, one might assume that the calculated t(∆) approximately coincide with the finite-slab Fresnel coefficients for a smooth material [63]. This produces an alternative way to extrapolate the index, which we find yields quantitatively similar results as what we present below. In Fig. 3, we plot our numerical results for the real and imaginary parts of n(∆), as a function of the input field detuning ∆, and for various densities. For comparison, we also plot the index as predicted by the MB equations, which starts to appreciably deviate from the full numerical results for dimensionless densities η 0.1. Interestingly, for sufficiently high densities, we observe that the computed spectra collapse onto the same curve when plotted as a function of the re-scaled detuning ∆/η, as shown in the insets of Fig. 3, which include all plots in the range 2 η 3. The invariance of n(∆/η) for η 2 directly indicates that both the maximum real index and the attenuation per unit length acquire fixed values with increasing density, and that density only determines a linear broadening in the spectra. Notably, the maximum real index saturates to a "real-life" value of ∼ 1.7, in contrast to the indefinite growth predicted by both MB and LL.
We note that a number of experiments involving dense cold atomic clouds have observed both a saturation of the index [38,40,41] and the emergence of an anomalous broadening of the linewidth [17,33,34,38,39], including a linear scaling with density [40,41]. A maximum index of n ≈ 1.26 has also been observed in experiments involving dense, hot atomic vapours [42], which has been attributed to atomic collisions. However, while complex collision dynamics necessitate semi-phenomenological models [64], here, our mechanism for saturation is quite fundamental, and occurs even for perfectly identical, stationary atoms.

III. RENORMALIZATION GROUP APPROACH
Our RG theory is based upon the key intuition gained in the collective scattering of just two atoms, to build up an understanding of the many-atom problem in a hierarchical manner. To be specific, let us consider the problem of two identical atoms, whose distance is much smaller than a wavelength, ρ 12 ≡ k 0 |r 1 − r 2 | 1. Applying Eq. 6 and Eq. 8, we can calculate the imaginary part of the "index" of the two-atom system, as illustrated in Fig. 4-a. One can see that the characteristic two-atom spectrum (blue line) is not twice the response of a single, isolated atom (green dashed curve), but instead consists of two, well-separated peaks with different linewidths and shifted resonances.
To understand this behavior, we consider the normal modes of the two-atom system, as encoded in the eigenstates of the dimensionless matrix G, whose elements G 12 were introduced in Eq. 7. When ρ 12 1, G is dominated by its off-diagonal components G 12 = G 21 , and in particular, by the purely real 1/ρ 3 12 near-field term (which we denote by G near 12 ). Specifically, in spherical coordinates ρ ij ≡ ρ ij (cos θ, sin θ sin φ, sin θ cos φ), G near ij = 3(−1 + 3 cos 2 θ)/(4ρ 3 ij ). This describes the strong, coherent, near-field coupling between the two dipoles. This produces symmetric and antisymmetric eigenstates whose dimensionless normal mode frequencies (real parts of the eigenvalues) are shifted as ω ± ≈ ±G near 12 , and align with the resonant peaks seen in Fig. 4-a. The linewidths (given by the imaginary parts of the eigenvalues) also are modified as Γ + ≈ 2 and Γ − ≈ ρ 2 12 , which is simply the two-atom limit of the famous Dicke superradiance model [23]. The key insight is that due to the large splitting, the total response in Fig. 4-a appears as two well-separated Lorentzians. Although this arose from the strong interaction of identical atoms, such a spectrum would also arise from two, inhomogeneous and non-interacting atoms, which were simply assigned these resonance frequencies and linewidths to start. This concept is at the heart of the RG approach for the many-atom case.
We now discuss how strong, coherent 1/ρ 3 ij near-field interactions in a many-atom system can be treated, by successively replacing strongly interacting pairs by optically equivalent, non-interacting atoms. Here, we will focus on the main conceptual steps of our RG scheme, while full details can be found in the Appendices. We anticipate that the scheme involves atoms with different renormalized resonant frequencies ω i , which can either interact, or not, through the near-field coupling, depending on the previous RG steps. The normal modes of such a system are given by the eigenstates of the generalized N × N matrix M = diag(ω) −G, where the elementsG ij are defined as Here, diag(ω) is a diagonal matrix containing the individual resonance frequencies ω = (ω 1 , . . . , ω N ), while L ij = 1 or 0 dictates whether pair i, j is allowed to interact via the near field (L ij = 1 for all pairs at the beginning of the RG process). In three dimensions, the 1/ρ 3 scaling of the near-field interaction implies that if an atom has a particularly close-by nearest neighbor, this pair will interact much more strongly between themselves than with any other nearby atoms [43]. Suppose that atoms i, j (with L ij = 1) are identified as the most strongly interacting pair, by a prescription given below. Then, we can re-write M as M = M pair + (M − M pair ), where the only non-zero elements of M pair involve atoms i, j. This effective 2 × 2 matrix reads where ω ij = (ω i + ω j )/2 and δω ij = (ω i − ω j )/2, and where we have included the coherent near-field interaction in M pair . The remaining far-field interactions between atoms i and j, as well as all the other atoms, are included in  To compare, we also plot twice the response of a single, isolated atom (green dashed line). b) Pictorial representation of the RG scheme. At each step of the RG flow the nearby pairs (identified by dashed circles) that mostly strongly interact via their near fields are identified, and replaced with atoms with different resonance frequencies (indicated by different colors) in such a way to produce an equivalent optical response. At the end of the RG process (last panel) the overall system is equivalent to an inhomogeneously broadened ensemble of weakly interacting atoms. c) Comparison between the maximum real refractive index predicted by the full coupled dipole simulations of identical atoms (blue points), and the equivalent, inhomogeneously broadened ensemble predicted by RG (green). For each value of density, the maximum index is obtained by optimizing over detuning. For comparison, the MB and LL models both predict a maximum index given by the orange curve. The inset compares the rescaled spectra Re n(∆/η) of the RG (green) and full coupled dipole (blue) simulations, given the points at densities η 2. d) Rescaled probability distribution of effective, inhomogeneously broadened resonance frequencies P (ωeff/η) obtained from the application of the RG scheme. Given 9 different values of the density η (ranging from η ≈ 2.5 up to η ≈ 80), the distributions of effective resonance frequencies are plotted with a different color, according to the bar on the right. The exact values chosen for the curves are emphasized by dotted white lines in the color-bar. The curves at η ≈ 2.5 and η ≈ 3 are calculated using the cylindrical system described in Fig. 2 From the structure of M pair , we define the pairwise interaction parameter K ij = L ij |G near ij |/(|δω ij |+1). A large value of K ij (which requires L ij = 1) implies that the strong near-field interaction is able to strongly split the original resonances, accounting for possible differences in resonance frequencies of the pair. We thus identify the most strongly interacting pair as that with the largest value of K ij , as pictorially depicted in the first panel of Fig. 4-b. Diagonalization of M pair results in two, new interacting resonance frequencies ω ± = ω ij ± δω 2 ij + (G near ij ) 2 . We can then obtain an approximately equivalent system by replacing the two original resonance frequencies ω i,j with the new values ω ± (second panel of Fig. 4-b). While the resulting modes are in principle delocalized between atoms i, j, to facilitate the RG, we randomly assign ω + to either atom i or j, while ω − is then assigned to the other atom (see SI on the issue of replacing atoms i, j with two new atoms placed at the midpoint of the original locations). This new system is described by a renormalized interaction matrix M eff = diag(ω eff ) −G eff , where ω eff = (ω 1 , . . . , ω + , . . . , ω − , . . . ω N ) contains the two renormalized resonance frequencies, and whereG eff includes the new set of allowed near-field interactions L eff , which both forbid the renormalized pair from interacting again (i.e. L eff ij = 0) and prevent any backflow of the RG process (see Appendix A for more details). The RG process can be iteratively repeated by identifying, at each step, the most strongly interacting pairs, and ends once K ij ≤ K cut-off ∼ 1, i.e. when all strong near-field interactions have been removed. In the numerics presented here, we take a cutoff parameter of K cut-off = 1. Other choices result in minor quantitative corrections, while the overall conclusions remain the same. The final result, as suggested in the third panel of Fig. 4-b, is that the original, homogeneous system can be mapped to an optically equivalent system that is inhomogeneously broadened, with a smooth probability distribution of resonance frequencies P (ω eff ).
Before describing more the characteristics of P (ω eff ), we discuss a few subtle but important points associated with the RG process. First, we point out the historic work of [43], which used RG to understand the properties of permanent, static dipoles, which only experience a near-field 1/ρ 3 interaction. Given only a near-field interaction in three dimensions, the interaction of a dipole with its nearest neighbor is then indeed dominant. However, we have a qualitatively different system, of driven, radiating dipoles. Naively then, a similar argument considering the 1/ρ far field would suggest that atoms within a shell of radius ρ and ρ + dρ of one atom at the origin would contribute an interaction strength of ∼ ρdρ, such that the furthest atoms actually play the strong role. We argue that an RG process based on the near field is still the correct prescription, as the index should be a local property. Instead, the apparent "dominance" of the far field simply reflects the fact that the macroscopic geometry of an optical system (e.g., if it is shaped as a lens or prism) can still drastically alter the overall optical response. Separately, we note that although the problem of just two atoms ( Fig. 4-b) can be interpreted in terms of renormalized resonance frequencies and linewidths, in the many-atom case, we only renormalize the resonance frequencies. Physically, this is because the decay rates of collective manyatom eigenstates reflect the rate of radiation of energy into the far field. This global property, which depends on the interference of all atoms, then generally cannot be derived from considering the isolated properties of individual pairs alone. Mathematically, this can be seen as the imaginary part of G ij (ρ ij ), which dictates the collective linewidths, does not contain a 1/ρ 3 ij near-field component. To validate the RG approach, we can use Eq. 6 and Eq. 7 (with the near-field interactions of renormalized atoms suitably removed, see Appendix A) to calculate the maximum real index as a function of density η of the ensemble with renormalized resonance frequencies. This is plotted in Fig. 4-c (green), along with exact numerical simulations (blue) of Eq. 6 for the original system of identical atoms. These curves show good agreement for all densities, and in particular, reveal a maximum index of n ≈ 1.7 at high densities. For comparison, the maximum index of the MB and LL equations (orange) increase indefinitely with density. Furthermore, motivated by our previous observation that high-density spectra collapse onto the same curve when the detuning is rescaled by density (insets of Fig. 3 and Fig. 4-c), in Fig. 4-d, we plot the re-scaled probability distribution of effective resonance frequencies P (ω eff /η) predicted by RG. For all densities considered (2.5 η 80), we see that a single universal curve results, i.e. the amount of broadening grows directly with density. Based on this curve, we find that the number of near-resonant atoms per reduced cubic wavelength (λ 0 /2π) 3 = k −3 0 , within a range ±Γ 0 of the original atomic resonance frequency, is approximately ∼ 0.3. The limited number of near-resonant atoms for light to interact with, regardless of how high the physical density is, directly explains the saturation of the maximum achievable index. We note that obtaining P (ω eff ) by RG does not require solving the coupled equations of Eq. 7, and we can calculate this distribution for much higher densities up to η ∼ 80. Furthermore, this distribution does not depend on the specific geometry. In Fig. 4-d, the curves for η ≤ 3 are obtained by a cylindrical geometry (the highest densities that we can compare to full coupled dipole simulations, as in Fig. 4-c), while for η > 3 we use a spherical geometry.
Within the language of RG, the universal distribution P (ω eff /η) constitutes the (numerically obtained) fixed point, as the interaction parameter of a system flows toward K ij → 1. While it might be desirable to write down and analytically solve the RG flow equation for P (ω eff ), this appears quite challenging in our case. This is because K ij not only depends on the distance between atoms, but also their spatial orientation (as the near field is anisotropic) and the difference in resonance frequencies.
As mentioned earlier, it is rather inconvenient to derive key optical properties of a system, like index, by solving a set of equations (Eq. 7) as large as the number of particles. At the same time, the RG approach clearly shows why conventional models (such as MB and LL) to treat atoms as a smooth medium fail at high densities [12,15,18], since the properties depend highly on granularity and the strong interaction between single nearest neighbors. Interestingly, RG also provides a basis to develop a more accurate smooth medium model. In particular, after the system is mapped to an inhomogeneously broadened distribution, P (ω eff ), where near-field interactions are seen to be strongly reduced, one can finally apply a smooth medium approximation. Specifically, Eq. 4 can be readily generalized to an inhomogeneously broadened ensemble Substituting the distribution found in Fig. 4-d, at high densities η 1, this equation predicts a maximum index of n ≈ 1.8, in good agreement with full results.

IV. DISCUSSION
To summarize, we have shown that despite the large resonant scattering cross section of a single atom, a dense atomic medium does not exhibit an anomalously large optical response. Rather, strong near-field interactions between atomic pairs combined with spatial disorder results in an effective inhomogeneous broadening mechanism, which occurs even if the atoms are otherwise perfect, and yields a maximum index of n ≈ 1.7. The key role of atomic granularity in this process also illustrates why conventional smooth medium approximations fail to describe the near-resonant response.
While we have focused on the linear refractive index, we believe that our RG formalism is valid in general for disordered atomic media, and constitutes a versatile new tool to study multiple scattering. Within the linear regime, RG might be used to provide additional insight to the question of whether an Anderson localization transition exists in a 3D ensemble, and under what conditions [31,32,53,[65][66][67]. Furthermore, it would be interesting to explore the usage of RG toward the challenging problem of quantum and nonlinear scattering. In the dilute limit, perturbative diagrammatic approaches [54] have only recently been developed. We hypothesize that a diagrammatic theory can also be developed in the dense, strong scattering regime, where strong interactions between nearby pairs are first nonperturbatively summed, while remaining interactions can be treated perturbatively.
Our results could also have interesting implications for quantum technologies based on atomic ensembles. In particular, the total optical depth of system, given by the product of the imaginary part of the index and system length, D ∼ (Im n)k 0 L, is a fundamental resource [68][69][70], with its magnitude establishing fundamental error bounds for most applications. As the imaginary part of the index also saturates with increasing density, this could place minimum size constraints on systems in order to achieve a given fidelity. Likewise, constraints on the maximum density could arise due to the induced inhomogeneous broadening, which typically constitutes an undesirable dephasing mechanism.
Finally, it would be interesting to understand more fully how the optical properties of a dilute atomic medium eventually transform into the low refractive index of actual optical materials, as the density is increased. Specifically, for a disordered ensemble, we have seen that the maximum index already saturates, at densities that are approximately six orders of magnitude before the onset of chemical processes. We hypothesize that the onset of chemistry, and the phase transition toward a real material, does not qualitatively alter the optical response, provided that the system remains disordered and the electrons tightly bound. Separately, it would be interesting to explore the same questions and transition for spatially ordered atomic systems, where RG breaks down and one expects very different qualitative behavior, due to the possibility of strong constructive and destructive interference in light scattering. To understand the transition to real materials, one must develop a theory that combines quantum chemistry and multiple scattering, which should be a rich avenue for future research.
At each step of the RG flow, we evaluate the list of couplings K ij = L ij |G near ij |/(|δω ij |+1) (where δω ij = (ω i − ω j )/2), ordering them from the largest to smallest in amplitude. Nominally, we should select the most strongly interacting pair and renormalize the pair properties, but the computational cost of this approach would be unfeasible for large atom number. Due to this reason, we start from the most strongly interacting pair (say, i, j), select it, and remove from the list all other pairs containing one of those atoms (e.g. i, k or j, k). We then proceed iteratively, until we select N RG most strongly interacting pairs. We choose N RG to be a small fraction of the total atom number N (approximately ∼ 2.5%), since the maximum number of possible disjoint pairs scales as N/2. Nevertheless, we have checked that the results are insensitive to different choices.
Given each pair (i, j) of the selected set, we diagonalize M pair , and define its eigenvalues as the new effective resonances ω ± = ω ij ± δω 2 ij + (G near ij ) 2 , where ω ij = (ω i + ω j )/2. We then substitute the initial frequencies (ω i , ω j ) with the new two effective resonances in ω, the order of the labels being chosen randomly.
We need to impose that the pair does not interact anymore through the near field, meaning that we must replace L old ij = 1 with L new ij = 0. At the same time, at any given stage of the RG flow, the resonance frequencies of any pair of effective atoms i and j might have been derived from a set of previous RG steps involving a set of atoms with indices {I } and {J }, respectively. If the sets {I } and {J } have some non-zero intersection, then atoms i and j must be omitted from a subsequent frequency renormalization step. Not doing this would violate the principle of RG, that we are integrating or "freezing" out the degrees of freedom with the strongest interactions. Numerically, we efficiently enforce this by replacing L new ik = L new jk = L new ki = L new kj = L old ik L old jk , ∀k, anytime a pair (i, j) is renormalized. Since L has (at any step) zero-valued diagonal elements, this directly ensures that L new ij = 0. After all atoms of the step have been renormalized, we re-evaluate the new set of K parameters, and repeat the scheme. When all pairs exhibit K ≤ K cut-off = 1, we stop the RG flow, obtaining an ensemble of N inhomogeneously  broadened atoms. Given a fixed value of the density η, we repeat this process for ≈ 100 different spatial configurations, in order to build up the final distribution P (ω eff ).
We extract the optical properties from the renormalized ensemble by applying Eq. 7 of the main text, modified in order to account for the the new N × N matrix M emerging from the RG scheme. This reads