---
abstract: '*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.*'
author:
- Francesco Andreoli
- 'Michael J. Gullans'
- 'Alexander A. High'
- Antoine Browaeys
- 'Darrick E. Chang'
bibliography:
- 'main.bib'
title: The maximum refractive index of an atomic medium
---

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\sim \mathcal O(1)$. While we typically utilize materials far from their natural electronic resonances, this even holds true close to resonance [@Rustgi1961OpticalEv; @Philipp1963OpticalSemiconductors; @Walker1964UltravioletDiamond; @Lamy1977OpticalUltraviolet; @Aspnes1983DielectricEV; @Warren1984OpticalMicrowave; @Papadopoulos1991OpticalDiamond; @Kitamura2007OpticalTemperature]. 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 $\sigma_{\text{sc}} \sim\lambda_0^2$ for photons resonant with an atomic transition of wavelength $\lambda_0$ ([Fig. \[fig:scatt\_cross\_sec\]]{}-a), which far exceeds both the physical size of the atom or the typical lattice constant of a solid ($\lambda_0\sim 1\upmu \text{m}$ for a typical optical transition, compared to the Bohr radius $a_0\sim 0.1 \text{nm}$).

Specifically, in standard theories [@Jackson1998ClassicalElectrodynamics; @Grynberg2010IntroductionOptics], the macroscopic index of an atomic medium ([Fig. \[fig:scatt\_cross\_sec\]]{}-b) is constructed from the product of the single-atom polarizability and the atomic density, and around resonance its value $n\sim \sqrt{N\lambda_0^3/V}$ extrapolates to a maximum of $ \sim 10^5$ at solid densities ([Fig. \[fig:scatt\_cross\_sec\]]{}-c). It is well-known that this argument neglects multiple scattering of light and photon-mediated dipole-dipole interactions [@Fleischhauer1999RadiativeFormula; @Javanainen2016LightOptics], and substantial work has been devoted to explore their effects on various optical phenomena, such as collective resonance shifts [@Manassah2012CooperativeShift; @Keaveney2012CooperativeThickness; @Javanainen2014ShiftsSample; @Bromley2016CollectiveMedium; @Jenkins2016OpticalGases; @Dobbertin2020CollectiveNanocavities], cooperative scattering properties [@Schilder2016PolaritonicAtoms; @Schilder2017HomogenizationScatterers; @Schilder2020Near-ResonantAtoms], emergence of sub- and super-radiance [@Dicke1954CoherenceProcesses; @Gross1982Superradiance:Emission; @Roof2016ObservationAtoms; @Araujo2016SuperradianceRegime; @Asenjo-Garcia2017ExponentialArrays; @He2019GeometricEmission], realization of atomic mirrors [@Shahmoon2017CooperativeArrays; @Bettles2016EnhancedArray; @Rui2020ALayer], and Anderson localization of light [@Skipetrov2014AbsenceScatterers; @Skipetrov2019SearchField]. In particular, this includes theoretical and experimental evidence that the optical response of dense gases can be much smaller than standard predictions [@Pellegrino2014ObservationEnsemble; @Schilder2016PolaritonicAtoms; @Jennewein2018CoherentTheory; @Javanainen2016LightOptics; @Guerin2017LightEffects] or even reach limiting values [@Chomaz2012AbsorptionAnalysis; @Zhu2016LightMedia; @Jenkins2016CollectiveExperiment; @Jennewein2016PropagationAtoms; @Corman2017TransmissionAtoms; @Jennewein2016CoherentCloud; @Keaveney2012MaximalNanolayer]. However, an underlying physical explanation is still missing, and our goal here is to understand better the mechanisms that might limit the index even when operating close to resonance.

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 $\sim 25000$ atoms, about an order of magnitude larger than comparable works [@Chomaz2012AbsorptionAnalysis; @Javanainen2014ShiftsSample; @Javanainen2016LightOptics; @Zhu2016LightMedia; @Schilder2016PolaritonicAtoms; @Schilder2017HomogenizationScatterers; @Schilder2020Near-ResonantAtoms; @Jennewein2016PropagationAtoms; @Jennewein2018CoherentTheory; @Corman2017TransmissionAtoms; @Jenkins2016CollectiveExperiment; @Dobbertin2020CollectiveNanocavities]), we find that the maximum index does not indefinitely grow with density, and saturates to a maximum value of $n\approx 1.7$, when the typical distance between atoms becomes smaller than the length scale associated with the resonant cross section, i.e. $d<\lambda_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 [@Levitov1990DelocalizationInteraction; @Fisher1994RandomChains; @Damle2000DynamicsChains; @Motrunich2000Infinite-randomnessPoints; @Refael2004EntanglementDimension; @Igloi2005StrongSystems; @Vosk2013Many-bodyPoint; @Refael2013StrongTransition]. In the context of our particular problem, the combination of strong near-field ($\sim 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 $\sim 1$ near-resonant atom per reduced cubic wavelength $\lambda_0^3/(2\pi)^3$, regardless of the physical atomic density, thus limiting the optical response ([Fig. \[fig:scatt\_cross\_sec\]]{}-d).

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 [@Chomaz2012AbsorptionAnalysis; @Lagendijk1996ResonantLight; @Skipetrov2014AbsenceScatterers; @Schilder2017HomogenizationScatterers; @Javanainen2016LightOptics; @Fayard2015IntensityPatterns; @Guerin2017LightEffects; @Corman2017TransmissionAtoms; @Zhu2016LightMedia; @Jennewein2018CoherentTheory; @Cottier2019MicroscopicLight], including in the nonlinear and quantum regimes [@Binninger2019NonlinearCloud].

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\approx 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.

![**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)\gg \lambda_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({\mathbf{r}},\omega_0)=|{\mathbf{E}}({\mathbf{r}},\omega_0)|^2$ (with the value indicated in the colorbar) over a hemispherical surface far from the ensemble (the radius of this hemisphere is $35\lambda_0$), given an input resonant Gaussian beam. The intensity is calculated for a single, random atomic configuration. The system parameters used are: beam waist $w_0=3\lambda_0$, cylinder radius $l_{\text{cyl.}}=7 \lambda_0$ and thickness $d=2\lambda_0$. []{data-label="fig:system_representation"}](FINAL_scatt_cross_sec.png){width="\columnwidth"}

![**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)\gg \lambda_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({\mathbf{r}},\omega_0)=|{\mathbf{E}}({\mathbf{r}},\omega_0)|^2$ (with the value indicated in the colorbar) over a hemispherical surface far from the ensemble (the radius of this hemisphere is $35\lambda_0$), given an input resonant Gaussian beam. The intensity is calculated for a single, random atomic configuration. The system parameters used are: beam waist $w_0=3\lambda_0$, cylinder radius $l_{\text{cyl.}}=7 \lambda_0$ and thickness $d=2\lambda_0$. []{data-label="fig:system_representation"}](FINAL_SM_system_and_speckles.pdf){width="\columnwidth"}

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 ${\left|g\right\rangle}$, ${\left|e\right\rangle}$, with frequency difference $\omega_0$ and associated wavelength $\lambda_0$, and which have an optically allowed transition with a dipole matrix element along a fixed axis (say ${\hat{\mathbf{x}}}$), as depicted in [Fig. \[fig:scatt\_cross\_sec\]]{}-a. The excited states of the atoms decay purely radiatively, with a rate of $\Gamma_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(\omega)$, we consider that the atoms are driven by a monochromatic, linearly-polarized input beam ${\mathbf{E}}_{\text{in}}({\mathbf{r}},\omega)=E_{\text{in}}({\mathbf{r}},\omega){\hat{\mathbf{x}}}$, whose polarization aligns with the polarizability axis of the atoms. Each atom $j$ acquires a dipole moment ${\mathbf{d}}_j(\omega)=d_j(\omega){\hat{\mathbf{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 [@Novotny2009PrinciplesNano-optics] [ $$\begin{array}{c}

\label{eq:FTMS_output_field}
{\mathbf{E}}({\mathbf{r}}, \omega)= {\mathbf{E}}_{\text{in}}({\mathbf{r}}, \omega) +\mu_0 \omega^2{
\displaystyle \sum 
}_{j=1}^N {{\Bar{\Bar{\mathbf{G}}}}}({\mathbf{r}}, {\mathbf{r}}_j, \omega)\cdot {\mathbf{d}}_j(\omega).

\end{array}$$ ]{} Here, the dyadic Green’s tensor ${{\Bar{\Bar{\mathbf{G}}}}}({\mathbf{r}}, {\mathbf{r}}_j, \omega)$ encodes the field at position ${\mathbf{r}}$, produced by an oscillating dipole at ${\mathbf{r}}_j$, and in vacuum is given by [@Novotny2009PrinciplesNano-optics] [ $$\begin{array}{c}

\label{eq:FTMS_Green_tensor}
{{\Bar{\Bar{\mathbf{G}}}}}({\mathbf{r}},{\mathbf{r}}',\omega) =
k \dfrac{e^{i \rho}}{4\pi }\left[\left(\dfrac{1}{\rho}+\dfrac{i}{\rho^2}-\dfrac{1}{\rho^3}\right)\mathbb I
+\left(-\dfrac{1}{\rho}-\dfrac{3i}{\rho^2}+\dfrac{3}{\rho^3}\right)\dfrac{\boldsymbol{\rho} \otimes \boldsymbol{\rho}  }{\rho^2}\right],

\end{array}$$ ]{} with $\rho \equiv |\boldsymbol\rho|\equiv k| ({\mathbf{r}} - {\mathbf{r}}')|$ and $k=\omega/c$. Note that ${{\Bar{\Bar{\mathbf{G}}}}}({\mathbf{r}},{\mathbf{r}}',\omega)$ contains both non-radiative, near-field ($\sim 1/\rho^3$) and radiative, far-field ($\sim 1/\rho$) terms.

Then, the induced dipole moment of atom $i$ is given by [ $$\begin{array}{c}

\label{eq:FTMS_coupled_dipoles}
 d_i(\omega)=\alpha_0(\omega) \epsilon_0\left[ E_{\text{in}}({\mathbf{r}}_i,\omega)+\mu_0 \omega^2 {
\displaystyle \sum 
}_{j\neq i}^{N-1} {\hat{\mathbf{x}}} \cdot {{\Bar{\Bar{\mathbf{G}}}}}({\mathbf{r}}_i,{\mathbf{r}}_j, \omega)\cdot {\hat{\mathbf{x}}}\;  d_j(\omega)\right],
\end{array}$$ ]{} where the parameter $\alpha_0(\omega)$ defines the polarizability of a single dipole. Although [Eq. \[eq:FTMS\_output\_field\]]{} and [Eq. \[eq:FTMS\_coupled\_dipoles\]]{} can describe any system of linearly-polarizable point-like dipoles [@Novotny2009PrinciplesNano-optics], e.g. dielectric nano-particles [@GarciaDeAbajo2007CollectiveMatter], in our case we focus on the response of non-absorbing, purely radiative atoms, whose resonant cross section $\sigma_{\text{sc}}=3\lambda_0^2/(2\pi)$ is the maximum set by the unitarity limit [@DeVries1998PointWaves]. In this context, the atomic polarizability reads $\alpha_0(\omega)=-3\pi /[(\Delta + i/2)k_0^3]$, where $k_0=2\pi/\lambda_0$ denotes the resonant wavevector, while $\Delta\equiv(\omega - \omega_0)/\Gamma_0$ represents the dimensionless detuning between the input beam frequency $\omega$ and the atomic resonance $\omega_0$.

While [Eq. \[eq:FTMS\_output\_field\]]{} and [Eq. \[eq:FTMS\_coupled\_dipoles\]]{} 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 [@Jackson1998ClassicalElectrodynamics] or equivalently the Maxwell-Bloch (MB) equations [@Grynberg2010IntroductionOptics], where the discreteness of atoms is replaced by a smooth medium of density $N/V$ ([Fig. \[fig:scatt\_cross\_sec\]]{}-c). The resulting index depends on the product of density and single-atom polarizability, [ $$\begin{array}{c}

\label{eq:FTMS_index_prediction_MB}
n_{\text{MB}}(\Delta)=\sqrt{1+\dfrac{N}{V}\alpha_0(\omega)}= \sqrt{1+ \dfrac{3\pi\eta}{-\Delta - i/2}},

\end{array}$$ ]{} where we defined the dimensionless density $\eta \equiv N/(Vk_0^3)$. Notably, for an optimum detuning, the maximum real part of the index scales as $\sim\sqrt{\eta}$.

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 ${\mathbf{r}}_0$, the model approximates the neighbouring atoms as a smooth dielectric medium with a small spherical exclusion around ${\mathbf{r}}_0$ [@Jackson1998ClassicalElectrodynamics]. The resulting local field correction gives an index that satisfies the equation $(n_{\text{LL}}^2-1)/(n_{\text{LL}}^2+2) = (N/V)\alpha_0(\omega)/3$ [@Novotny2009PrinciplesNano-optics]. Plugging in the atomic polarizability, one readily finds that [ $$\begin{array}{c}

\label{eq:FTMS_index_prediction_LL}
n_{\text{LL}}(\Delta)=n_{\text{MB}}(\Delta+\pi\eta).
\end{array}$$ ]{} Importantly, while the spectrum is shifted, the LL model still produces a maximum index that grows like $\sim\sqrt{\eta}$.

![ **Frequency-dependent refractive index for different atomic densities.** The solid lines portray the imaginary (subfigure $a$) and real (subfigure $b$) part of the refractive index versus dimensionless detuning $\Delta$, obtained through [Eq. \[eq:CDS\_input-output\_eq\_single\_mode\]]{}, while the dotted lines show the MB predictions. The colors denote different atomic densities (colorbar on right), with the specific values indicated by the dotted white lines. The refractive index is inferred by averaging the complex transmission coefficient $t(\Delta)$ over $\sim 10^3 - 10^4$ atomic configurations. Other system parameters are: thickness $d =0.4\lambda_0$, transverse radius $5 \leq l_{\text{cyl.}}/\lambda_0\leq 7$, beam waist $2.5\leq w_0/\lambda_0 \leq 3$. The insets show the curves at the 3 highest densities as a function of the rescaled detuning $\Delta/\eta$. []{data-label="fig:SM_spectrum"}](FINAL_SM_spectrum.pdf){width="\columnwidth"}

Coupled-dipole simulations
==========================

[Eq. \[eq:FTMS\_output\_field\]]{} and [Eq. \[eq:FTMS\_coupled\_dipoles\]]{} 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. \[fig:system\_representation\]]{}-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\lesssim \lambda_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 [@Shapiro1986LargeMedia; @Lagendijk1996ResonantLight; @Fayard2015IntensityPatterns; @Schilder2017HomogenizationScatterers] that for a fixed random spatial configuration, an input as in [Fig. \[fig:system\_representation\]]{} 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. \[fig:system\_representation\]]{}-b. 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. \[eq:FTMS\_output\_field\]]{} 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(\Delta)$ given by [@Chomaz2012AbsorptionAnalysis; @Manzoni2018OptimizationArrays] [ $$\begin{array}{c}

\label{eq:CDS_input-output_eq_single_mode}
t(\Delta)=1+\dfrac{3 i}{(w_0k_0)^2 }{
\displaystyle \sum 
}_{j=1}^N \dfrac{  E^*_{\text{in}}({\mathbf{r}}_j,\omega_0)}{E_0}  c_j(\Delta),
\end{array}$$ ]{} where $E_0$ is the input field amplitude at the beam focus. Here, for convenience, we have defined re-scaled dipole amplitudes $c_j(\Delta) = d_j(\omega) k_0^3/(3\pi \epsilon_0 E_0) $, which satisfy the dimensionless coupled equations [ $$\begin{array}{c}

\label{eq:CDS_input-output_eq_single_mode2}
-\Delta c_i(\Delta)-\displaystyle
{
\displaystyle \sum 
}_{j=1}^{N}  G_{ij} c_j(\Delta)   = \; \dfrac{ E_{\text{in}}({\mathbf{r}}_i,\omega_0)}{E_0},

\end{array}$$ ]{} In these equations, we define $G_{ij}\equiv ( 3\pi/k_0) {\hat{\mathbf{x}}} \cdot {{\Bar{\Bar{\mathbf{G}}}}}({\mathbf{r}}_i,{\mathbf{r}}_j,\omega_0)\cdot {\hat{\mathbf{x}}}$ and $G_{jj}=i/2$, which coincides with the single-atom decay rate in units of $\Gamma_0$, while regularizing the divergent self-energy associated with the real part of ${{\Bar{\Bar{\mathbf{G}}}}}$. Note that, for simplicity, the Green’s function ${{\Bar{\Bar{\mathbf{G}}}}}({\mathbf{r}}_i,{\mathbf{r}}_j,\omega_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 ${{\Bar{\Bar{\mathbf{G}}}}}$ is an excellent approximation for near-resonant atoms, as the optical dispersion and delay of such a system is dominated by the atomic response itself rather than from the vacuum [@Chang2012CavityMirrors]. Similarly, we approximate the near-resonant input field as $E_{\text{in}}({\mathbf{r}}_i,\omega)\simeq E_{\text{in}}({\mathbf{r}}_i,\omega_0)$.

The expression in [Eq. \[eq:CDS\_input-output\_eq\_single\_mode\]]{} represents a useful closed-form definition of the transmission coefficient $t(\Delta)$, which avoids a numerically expensive point-by-point evaluation of the scattered field ${\mathbf{E}}({\mathbf{r}},\omega)$, as nominally prescribed by [Eq. \[eq:FTMS\_output\_field\]]{}. We can extrapolate the complex index of refraction $n(\Delta)$ from the relation [ $$\begin{array}{c}

\label{eq:CDS_refractive_index_def_from_t}
{\langle t(\Delta) \rangle}= \exp \left\{ i  \left[n\left(\Delta \right)-1\right] k_0 d\right\}.

\end{array}$$ ]{} where the averages are performed over $\sim 10^3 - 10^4$ sets of random positions, for each fixed density. Unlike in a smooth medium, we have that $|{\langle  t(\Delta) \rangle}|^2 \neq {\langle |t(\Delta)|^2 \rangle}$. 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 [@Davis1995Bose-EinsteinAtoms; @Anderson1995ObservationVapor]). In the Appendix B, we discuss the independence of the calculated index from the thickness $d$, which is implicitly assumed in [Eq. \[eq:CDS\_refractive\_index\_def\_from\_t\]]{}. Alternatively, one might assume that the calculated ${\langle  t(\Delta) \rangle}$ approximately coincide with the finite-slab Fresnel coefficients for a smooth material [@Born1999PrinciplesOptics]. This produces an alternative way to extrapolate the index, which we find yields quantitatively similar results as what we present below.

In [Fig. \[fig:SM\_spectrum\]]{}, we plot our numerical results for the real and imaginary parts of $n(\Delta)$, as a function of the input field detuning $\Delta$, 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 $\eta \gtrsim 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 $\Delta/\eta$, as shown in the insets of [Fig. \[fig:SM\_spectrum\]]{}, which include all plots in the range $ 2 \lesssim \eta\lesssim 3$. The invariance of $n(\Delta/\eta)$ for $\eta \gtrsim 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 $\sim 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 [@Jenkins2016CollectiveExperiment; @Jennewein2016CoherentCloud; @Corman2017TransmissionAtoms] and the emergence of an anomalous broadening of the linewidth [@Pellegrino2014ObservationEnsemble; @Jenkins2016CollectiveExperiment; @Jenkins2016OpticalGases; @Jennewein2016PropagationAtoms; @Jennewein2018CoherentTheory], including a linear scaling with density [@Jennewein2016CoherentCloud; @Corman2017TransmissionAtoms]. A maximum index of $n\approx 1.26$ has also been observed in experiments involving dense, hot atomic vapours [@Keaveney2012MaximalNanolayer], which has been attributed to atomic collisions. However, while complex collision dynamics necessitate semi-phenomenological models [@Allard1982TheLines], here, our mechanism for saturation is quite fundamental, and occurs even for perfectly identical, stationary atoms.

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, $\rho_{12} \equiv k_0 |{\mathbf{r}}_1 - {\mathbf{r}}_2|\ll 1$. Applying [Eq. \[eq:CDS\_input-output\_eq\_single\_mode\]]{} and [Eq. \[eq:CDS\_refractive\_index\_def\_from\_t\]]{}, we can calculate the imaginary part of the “index" of the two-atom system, as illustrated in [Fig. \[fig:RG\_plots\]]{}-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. \[eq:CDS\_input-output\_eq\_single\_mode2\]]{}. When $\rho_{12}\ll 1$, $G$ is dominated by its off-diagonal components $G_{12}=G_{21}$, and in particular, by the purely real $1/\rho_{12}^3$ near-field term (which we denote by $G_{12}^{\text{near}}$). Specifically, in spherical coordinates $\boldsymbol \rho_{ij} \equiv \rho_{ij} (\cos \theta, \sin \theta \sin \phi, \sin \theta \cos \phi)$, $G_{ij}^{\text{near}}= 3(-1+3\cos^2 \theta)/(4\rho_{ij}^3) $. This describes the strong, coherent, near-field coupling between the two dipoles. This produces symmetric and anti-symmetric eigenstates whose dimensionless normal mode frequencies (real parts of the eigenvalues) are shifted as $\omega_{\pm}\approx \pm G_{12}^{\text{near}}$, and align with the resonant peaks seen in [Fig. \[fig:RG\_plots\]]{}-a. The linewidths (given by the imaginary parts of the eigenvalues) also are modified as $\Gamma_{+}\approx 2$ and $\Gamma_{-}\approx \rho_{12}^2$, which is simply the two-atom limit of the famous Dicke superradiance model [@Gross1982Superradiance:Emission]. The key insight is that due to the large splitting, the total response in [Fig. \[fig:RG\_plots\]]{}-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/\rho_{ij}^3$ 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 $\omega_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 \times N$ matrix $\mathcal M= \text{diag}( \boldsymbol \omega )-\tilde G$, where the elements $\tilde G_{ij}$ are defined as $\tilde G_{ij}= \mathcal L_{ij} G_{ij}^{\text{near}} + (G_{ij}-G_{ij}^{\text{near}})$. Here, $\text{diag}( \boldsymbol \omega )$ is a diagonal matrix containing the individual resonance frequencies $\boldsymbol \omega =(\omega_1,\dots,\omega_N)$, while $\mathcal L_{ij}=1$ or $0$ dictates whether pair $i,j$ is allowed to interact via the near field ($\mathcal L_{ij}=1$ for all pairs at the beginning of the RG process). In three dimensions, the $1/\rho^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 [@Levitov1990DelocalizationInteraction]. Suppose that atoms $i,j$ (with $\mathcal L_{ij}=1$) are identified as the most strongly interacting pair, by a prescription given below. Then, we can re-write $\mathcal M$ as $\mathcal M=\mathcal M_{\text{pair}}+(\mathcal M-\mathcal M_{\text{pair}})$, where the only non-zero elements of $\mathcal M_{\text{pair}}$ involve atoms $i,j$. This effective $2\times 2$ matrix reads [ $$\begin{array}{c}
\mathcal M_{\text{pair}} ={\langle \omega \rangle}_{ij} \mathbb{I}
+
\left(\begin{array}{cc}
\delta\omega_{ij} & - G^{\text{near}}_{ij}\\\\
-  G^{\text{near}}_{ij} & -\delta\omega_{ij}
\end{array}
\right),
\end{array}$$ ]{} where ${\langle \omega \rangle}_{ij}= (\omega_i +\omega_j)/2$ and $\delta\omega_{ij}= (\omega_i -\omega_j)/2$, and where we have included the coherent near-field interaction in $\mathcal M_{\text{pair}}$. The remaining far-field interactions between atoms $i$ and $j$, as well as all the other atoms, are included in $(\mathcal M-\mathcal M_{\text{pair}})$. The large near-field interaction motivates diagonalizing $\mathcal M_{\text{pair}}$ first, while treating $(\mathcal M-\mathcal M_{\text{pair}})$ as a perturbation.

![**Renormalization group analysis.** a) Representative optical response of two identical atoms separated by a distance $\rho_{12} \ll 1$. Here, we plot the absorption spectrum (blue curve), which consists of two well-separated Lorentzians. The positions of the resonances are given by $\sim \pm G_{12}^{\text{near}}$, where $G_{12}^{\text{near}}\propto 1/\rho_{12}^3$ is the near-field component of the Green’s function. 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 ${
\text{Re}\;n
}(\Delta/\eta)$ of the RG (green) and full coupled dipole (blue) simulations, given the points at densities $\eta \gtrsim 2$. d) Rescaled probability distribution of effective, inhomogeneously broadened resonance frequencies $P(\omega_{\text{eff}}/\eta)$ obtained from the application of the RG scheme. Given 9 different values of the density $\eta$ (ranging from $\eta \approx 2.5$ up to $\eta \approx 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 $\eta \approx 2.5$ and $\eta \approx 3$ are calculated using the cylindrical system described in [Fig. \[fig:system\_representation\]]{}-a (with thickness $d =0.4\lambda_0$ and transverse radius $l_{\text{cyl.}}= 5\lambda_0$), while the others are evaluated for a sphere of radius $r=0.55\lambda_0$. []{data-label="fig:RG_plots"}](FINAL_RG_plots.pdf){width="\columnwidth"}

From the structure of $\mathcal M_{\text{pair}}$, we define the pairwise interaction parameter $\mathcal K_{ij}=\mathcal L_{ij}|G^{\text{near}}_{ij}|/(|\delta \omega_{ij}|+1)$. A large value of $\mathcal K_{ij}$ (which requires $\mathcal 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 $\mathcal K_{ij}$, as pictorially depicted in the first panel of [Fig. \[fig:RG\_plots\]]{}-b. Diagonalization of $\mathcal M_{\text{pair}}$ results in two, new interacting resonance frequencies $\omega_{\pm}={\langle \omega \rangle}_{ij}\pm \sqrt{\delta\omega_{ij}^2 + (G_{ij}^{\text{near}})^2}$. We can then obtain an approximately equivalent system by replacing the two original resonance frequencies $\omega_{i,j}$ with the new values $\omega_{\pm}$ (second panel of [Fig. \[fig:RG\_plots\]]{}-b). While the resulting modes are in principle delocalized between atoms $i,j$, to facilitate the RG, we randomly assign $\omega_{+}$ to either atom $i$ or $j$, while $\omega_{-}$ 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 $\mathcal M_{\text{eff}}=\text{diag}(\boldsymbol{\omega}_{\text{eff}})-\tilde G_{\text{eff}}$, where $\boldsymbol{\omega}_{\text{eff}}=(\omega_1,\dots,\omega_{+}, \dots, \omega_{-},\dots \omega_N)$ contains the two renormalized resonance frequencies, and where $\tilde G_{\text{eff}}$ includes the new set of allowed near-field interactions $\mathcal L_{\text{eff}}$, which both forbid the renormalized pair from interacting again (i.e. $\mathcal L^{\text{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 $\mathcal K_{ij} \leq \mathcal K_{\text{cut-off}}\sim 1$, i.e. when all strong near-field interactions have been removed. In the numerics presented here, we take a cutoff parameter of $\mathcal K_{\text{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. \[fig:RG\_plots\]]{}-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(\omega_{\text{eff}})$.

Before describing more the characteristics of $P(\omega_{\text{eff}})$, we discuss a few subtle but important points associated with the RG process. First, we point out the historic work of [@Levitov1990DelocalizationInteraction], which used RG to understand the properties of permanent, *static* dipoles, which only experience a near-field $1/\rho^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/\rho$ far field would suggest that atoms within a shell of radius $\rho$ and $\rho+d\rho$ of one atom at the origin would contribute an interaction strength of $\sim \rho d\rho$, 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. \[fig:RG\_plots\]]{}-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 many-atom 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}(\rho_{ij})$, which dictates the collective linewidths, does not contain a $1/\rho_{ij}^3$ near-field component.

To validate the RG approach, we can use [Eq. \[eq:CDS\_input-output\_eq\_single\_mode\]]{} and [Eq. \[eq:CDS\_input-output\_eq\_single\_mode2\]]{} (with the near-field interactions of renormalized atoms suitably removed, see Appendix A) to calculate the maximum real index as a function of density $\eta$ of the ensemble with renormalized resonance frequencies. This is plotted in [Fig. \[fig:RG\_plots\]]{}-c (green), along with exact numerical simulations (blue) of [Eq. \[eq:CDS\_input-output\_eq\_single\_mode\]]{} for the original system of identical atoms. These curves show good agreement for all densities, and in particular, reveal a maximum index of $n\approx 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. \[fig:SM\_spectrum\]]{} and [Fig. \[fig:RG\_plots\]]{}-c), in [Fig. \[fig:RG\_plots\]]{}-d, we plot the re-scaled probability distribution of effective resonance frequencies $P(\omega_{\text{eff}}/\eta)$ predicted by RG. For all densities considered ($2.5 \lesssim \eta \lesssim 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 $(\lambda_0/2\pi)^3=k_0^{-3}$, within a range $\pm\Gamma_0$ of the original atomic resonance frequency, is approximately $\sim 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(\omega_{\text{eff}})$ by RG does not require solving the coupled equations of [Eq. \[eq:CDS\_input-output\_eq\_single\_mode2\]]{}, and we can calculate this distribution for much higher densities up to $\eta\sim 80$. Furthermore, this distribution does not depend on the specific geometry. In [Fig. \[fig:RG\_plots\]]{}-d, the curves for $\eta\leq 3$ are obtained by a cylindrical geometry (the highest densities that we can compare to full coupled dipole simulations, as in [Fig. \[fig:RG\_plots\]]{}-c), while for $\eta>3$ we use a spherical geometry.

Within the language of RG, the universal distribution $P(\omega_{\text{eff}}/\eta)$ constitutes the (numerically obtained) fixed point, as the interaction parameter of a system flows toward $\mathcal K_{ij}\rightarrow 1$. While it might be desirable to write down and analytically solve the RG flow equation for $P(\omega_{\text{eff}})$, this appears quite challenging in our case. This is because $\mathcal 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. \[eq:CDS\_input-output\_eq\_single\_mode2\]]{}) 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 [@Javanainen2014ShiftsSample; @Javanainen2016LightOptics; @Dobbertin2020CollectiveNanocavities], 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(\omega_{\text{eff}})$, where near-field interactions are seen to be strongly reduced, one can finally apply a smooth medium approximation. Specifically, [Eq. \[eq:FTMS\_index\_prediction\_MB\]]{} can be readily generalized to an inhomogeneously broadened ensemble [ $$\begin{array}{c}
n(\Delta)=\sqrt{1+3\pi \eta {
\displaystyle \int
}\dfrac{ P(\omega_{\text{eff}})}{-\Delta + \omega_{\text{eff} }- i/2} d\omega_{\text{eff}}}.
\end{array}$$ ]{} Substituting the distribution found in [Fig. \[fig:RG\_plots\]]{}-d, at high densities $\eta \gg 1$, this equation predicts a maximum index of $n\approx 1.8$, in good agreement with full results.

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\approx 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 [@Skipetrov2014AbsenceScatterers; @Bellando2014CooperativeHamiltonian; @Skipetrov2016RedLocalization; @Skipetrov2018Ioffe-RegelScatterers; @Cottier2019MicroscopicLight; @Skipetrov2019SearchField]. 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 [@Binninger2019NonlinearCloud] 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 non-perturbatively 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\sim ({
\text{Im}\;n
})k_0 L$, is a fundamental resource [@Gorshkov2007UniversalMedia; @Hammerer2010QuantumEnsembles; @Gorshkov2011Photon-PhotonBlockade], 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.

Acknowledgements {#acknowledgements .unnumbered}
================

We acknowledge S. Grava, N. Fayard, J.-J. Greffet, and H.J. Kimble for stimulating discussions. F.A. acknowledges support from the ICFOstepstone - PhD Programme funded by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 713729. D.E.C. acknowledges support from MINECO Severo Ochoa Grant SEV-2015-0522, CERCA Programme / Generalitat de Catalunya, Fundació Privada Cellex, Fundació Mir-Puig, ERC Starting Grant FOQAL, Plan Nacional Grant ALIQS (funded by MCIU, AEI, and FEDER), Fundación Ramón Areces Project CODEC, and QuantumCAT (funded within the framework of the ERDF Operational Program of Catalonia).

![**Independence of the refractive index from the thickness of the ensemble.** Given the physical system of Fig. 2-a of the main text (with $w_0=2.5\lambda_0$, $l_{\text{cyl.}}=5\lambda_0$), we compare the resonant ($\Delta = 0$) refractive index as a function of the density, for various ensemble thicknesses: $d=0.4\lambda_0$ (as in the main text, here in blue), $d=0.6\lambda_0$ (in green) and $d=0.8\lambda_0$ (in orange). Subfigures a) and b) illustrate the real and imaginary parts of the index, respectively. The insets show the full spectra $n(\Delta)$ at a fixed density $\eta \simeq 0.28$. All data are obtained by averaging ${\langle  t(\Delta) \rangle}$ over $>1000$ configurations. []{data-label="fig:SI_thick_indep"}](FINAL_SI_thick_indep.pdf){width="\columnwidth"}

Appendix A: Full description of the RG algorithm {#appendix-a-full-description-of-the-rg-algorithm .unnumbered}
------------------------------------------------

Here, we provide a full description of the RG algorithm. We assume that we have an ensemble of $N$ randomly positioned atoms. As shown in [Eq. \[eq:CDS\_input-output\_eq\_single\_mode2\]]{}, each pair of atoms interacts through the coupling $G_{ij}= ( 3\pi/k_0) {\hat{\mathbf{x}}} \cdot {{\Bar{\Bar{\mathbf{G}}}}}(\boldsymbol \rho_{ij},\omega_0)\cdot {\hat{\mathbf{x}}}$, where $\boldsymbol \rho_{ij} \equiv k_0 ({\mathbf{r}}_i-{\mathbf{r}}_j)$. The $1/\rho_{ij}^3$ near-field component of $G_{ij}$ reads [ $$\begin{array}{c}

G_{ij}^{\text{near}}= \dfrac{3}{4\rho_{ij}^3}(-1+3\cos^2 \theta) ,

\end{array}$$ ]{} where we have represented $\boldsymbol \rho_{ij} \equiv \rho_{ij} (\cos \theta, \sin \theta \sin \phi, \sin \theta \cos \phi)$ in spherical coordinates. This interaction is purely real, and describes a coherent interaction between dipoles.

Let us now consider a generic step of the RG flow, where the atomic ensemble is already composed of effective atoms characterized by different atomic resonances and a specific set of allowed near-field interactions. As discussed in the main text, this system is described by the $N\times N$ matrix $\mathcal M= \text{diag}( \boldsymbol \omega )-\tilde G$, where the elements $\tilde G_{ij}$ read $\tilde G_{ij}= \mathcal L_{ij} G_{ij}^{\text{near}} + (G_{ij}-G_{ij}^{\text{near}})$. Numerically, this matrix is initialized according to $\boldsymbol \omega ^{\text{init.}}=(0,\dots,0)$ and $\mathcal L_{ij}^{\text{init.}}=1-\delta_{ij}$, stating that all atoms are resonant at the frequency $\omega_0$ and cannot self-interact.

At each step of the RG flow, we evaluate the list of couplings $\mathcal K_{ij}=\mathcal L_{ij}|G^{\text{near}}_{ij}|/(|\delta \omega_{ij}|+1)$ (where $\delta\omega_{ij}= (\omega_i -\omega_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_{\text{RG}}$ most strongly interacting pairs. We choose $N_{\text{RG}}$ to be a small fraction of the total atom number $N$ (approximately $\sim 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 $\mathcal M_{\text{pair}} $, and define its eigenvalues as the new effective resonances $\omega_{\pm}={\langle \omega \rangle}_{ij}\pm \sqrt{\delta\omega_{ij}^2 + (G_{ij}^{\text{near}})^2}$, where ${\langle \omega \rangle}_{ij}= (\omega_i +\omega_j)/2$. We then substitute the initial frequencies ($\omega_i$, $\omega_j$) with the new two effective resonances in $\boldsymbol{\omega}$, 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 $\mathcal L_{ij}^{\text{old}}=1$ with $\mathcal L_{ij}^{\text{new}}=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 $\mathcal L_{ik}^{\text{new}}=\mathcal L_{jk}^{\text{new}}=\mathcal L_{ki}^{\text{new}}=\mathcal L_{kj}^{\text{new}}=\mathcal L_{ik}^{\text{old}}\mathcal L_{jk}^{\text{old}},$ $\forall k$, anytime a pair $(i,j)$ is renormalized. Since $\mathcal L$ has (at any step) zero-valued diagonal elements, this directly ensures that $\mathcal L_{ij}^{\text{new}}=0$.

After all atoms of the step have been renormalized, we re-evaluate the new set of $\mathcal K$ parameters, and repeat the scheme. When all pairs exhibit $\mathcal K\leq \mathcal K_{\text{cut-off}}=1$, we stop the RG flow, obtaining an ensemble of $N$ inhomogeneously broadened atoms. Given a fixed value of the density $\eta$, we repeat this process for $\approx 100$ different spatial configurations, in order to build up the final distribution $P(\omega_{\text{eff}})$.

We extract the optical properties from the renormalized ensemble by applying [Eq. \[eq:CDS\_input-output\_eq\_single\_mode2\]]{} of the main text, modified in order to account for the the new $N\times N$ matrix $\mathcal{M}$ emerging from the RG scheme. This reads

\[eq:Methods\_input-output\_eq\_single\_mode\] (-+ \_i) c\_i()-\_[j,i=1]{}\^[N]{} c\_j() = .

Appendix B: Linear behaviour of the index as a function of thickness. {#appendix-b-linear-behaviour-of-the-index-as-a-function-of-thickness. .unnumbered}
---------------------------------------------------------------------

Our operative definition of the complex index of refraction, as given by Eq. 7 of the main text, is [ $$\begin{array}{c}

{\langle t(\Delta) \rangle}= \exp \left\{ i  \left[n\left(\Delta \right)-1\right] k_0 d\right\}.

\end{array}$$ ]{} Since the refractive index is an intensive property by definition, it must not depend upon the thickness $d$ that we choose in our numerics. Here, we show that our operative definition satisfies this condition.

We consider the same physical system described in Fig. 2-a of the main text, with $w_0=2.5\lambda_0$, $l_{\text{cyl.}}=5\lambda_0$ and different values of the thickness. By applying Eq. 6 and Eq. 7 of the main text, we compute the resonant ($\Delta = 0$) refractive index for growing values of the density $\eta$, and we plot its real (imaginary) part in [Fig. \[fig:SI\_thick\_indep\]]{}-a(b). The simulated values of the thickness are: $d=0.4\lambda_0$ (as in the main text, here in blue), $d=0.6\lambda_0$ (in green) and $d=0.8\lambda_0$ (in orange). Moreover, for the point at $\eta\simeq 0.28$, we evaluate the full spectra $n(\Delta)$, as represented in the insets of the figure. All curves show the same behaviour, independently of $d$, both on resonance and varying the detuning.

Appendix C: Definition of effective positions in the RG scheme. {#appendix-c-definition-of-effective-positions-in-the-rg-scheme. .unnumbered}
---------------------------------------------------------------

In the main text (cf. Fig. 4b), we described how the optical response of a pair of atoms separated by a distance $\rho_{ij}\ll 1$ is characterized by two effective resonance frequencies, corresponding to the real parts of the eigenvalues of the two-atom system. The two collective modes are intrinsically delocalized in space (being formed by atoms with two different positions ${\mathbf{r}}_{i,\;j}$). As this delocalization is difficult to incorporate into the RG scheme, we instead attribute each of these two resonance frequencies to a new effective atom, with well-defined position.

In the main text, it was stated that the new effective atomic positions are assigned to those of the original pair, ${\mathbf{r}}_{i,\;j}$ (randomly between the two possible permutations). A more natural choice, given that the two renormalized atoms are non-interacting, might be to place them at the midpoint $({\mathbf{r}}_{i}+{\mathbf{r}}_{j})/2$ between the two original atoms, but here we discuss the problem with that approach.

Specifically, for a finite-size sample, the atoms closest to the perimeter of the sample will only renormalize with atoms that are closer to the interior. As illustrated in [Fig. \[fig:SI\_shrink\]]{}, this means that step by step, the shape of the cloud tends to shrink. This effectively distorts the ensemble and results in a higher density, and higher interaction strengths in the next step of RG.

![**RG scheme based upon re-positioning atoms.** Due to the finite size of the sample, if one defines the positions of the new effective atoms as being at the midpoint between the original pair, then, at each RG step, the cloud effectively shrinks, resulting in a distortion of the ensemble.[]{data-label="fig:SI_shrink"}](FINAL_SI_shrink.png){width="0.75\columnwidth"}
