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. 