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 