# LYRA III: The smallest Reionization survivors

Thales A. Gutcke Department of Astrophysical Sciences, Princeton University, 4
Ivy Lane, Princeton, NJ 08544, USA  Christoph Pfrommer Leibniz-Institute for
Astrophysics Potsdam (AIP), An der Sternwarte 16, 14482 Potsdam, Germany  Greg
L. Bryan Center for Computational Astrophysics, Flatiron Institute, 162 5th
Ave, New York, NY 10010, USA  Columbia University, Department of Astronomy,
550 W 120th St, New York, NY, 10025, USA  Rüdiger Pakmor Max-Planck-Institut
für Astrophysik, Karl-Schwarzschild-Str. 1, D-85748, Garching, Germany  Volker
Springel Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1,
D-85748, Garching, Germany  Thorsten Naab Max-Planck-Institut für Astrophysik,
Karl-Schwarzschild-Str. 1, D-85748, Garching, Germany

###### Abstract

The dividing line between galaxies that are quenched by reionization
(“relics”) and galaxies that survive reionization (i.e. continue forming
stars) is commonly discussed in terms of a halo mass threshold. We probe this
threshold in a physically more complete and accurate way than has been
possible to date, using five extremely high resolution
($M_{\mathrm{target}}=4~{}\,\mathrm{M}_{\odot}$) cosmological zoom-in
simulations of dwarf galaxies within the halo mass range $1-4\times
10^{9}~{}\,\mathrm{M}_{\odot}$. The employed LYRA simulation model features
resolved interstellar medium physics and individual, resolved supernova
explosions. In our results, we discover an interesting intermediate population
of dwarf galaxies close to the threshold mass but which are neither full
reionization relics nor full reionization survivors. These galaxies initially
quench at the time of reionization but merely remain quiescent for $\sim 500$
Myr. At $z\sim 5$ they recommence star formation in a synchronous way, and
remain star-forming until the present day. These results demonstrate that the
halo mass at $z=0$ is not a good indicator of survival close to the threshold.
While the star formation histories we find are diverse, we show that they are
directly related to the ability of a given halo to retain and cool gas.
Whereas the latter is most strongly dependent on the mass (or virial
temperature) of the host halo at the time of reionization, it also depends on
its growth history, the UV background (and its decrease at late times) and the
amount of metals retained within the halo.

††thanks: NASA Hubble Fellow

##  1 Introduction

Dwarf galaxies are observed to display a variety of star formation histories
(SFHs, Weisz et al. 2014). Some SFHs show fairly continuous activity to the
present day. Others are quenched at high redshift and then “rejuvenated” in
the last few gigayears before the present day. Lastly, there are “reionization
relics” that formed all their stars before the Universe was reionized ($z\sim
8-7$) and have remained “red and dead” since. To obtain SFHs from
observations, they can be reconstructed either from the color-magnitude
diagram (CMD) or the spectral energy distribution (SED). This has been done
for known dwarf galaxies in and around the Local Group, to a maximum distance
of $\sim 4~{}\rm{Mpc}$ (Weisz et al., 2014; McQuinn et al., 2010; Olsen et
al., 2021). Olsen et al. (2021) specifically find synchronized SF across their
sample, requiring a large scale environmental cause or a general cosmological
explanation. However, observations of such small systems are still only
possible very close to the Local Group. Disentangling the contribution of
reionization and environment is thus not easily possible for such faint
systems.

From a theoretical point of view, there should be a cut-off mass for galaxy
formation due to the heating effect of the photoionizing background (Rees,
1986; Efstathiou, 1992; Babul & Rees, 1992). This was simulated in detail by
various authors (e.g. Thoul & Weinberg, 1996; Gnedin, 2000; Benson et al.,
2002; Hoeft et al., 2006; Okamoto et al., 2008; Hambrick et al., 2009), who
find substantial gas mass loss and suppression of SF in low mass halos due the
heating effect by the ultra-violet background (UVB). Okamoto et al. (2008)
predict the characteristic, cut-off halo mass to be $M_{c}\sim 6.5\times
10^{9}~{}\,\mathrm{M}_{\odot}$ at $z=0$.

However, the characteristic mass is dependent on redshift and is lower at
higher redshift. Macciò et al. (2017) distinguish halos that remain dark or
manage to form stars and propose that the ones that remain dark do so because
their halos do not cross the characteristic mass at any time in their
evolution. Dijkstra et al. (2004) revisit the Thoul & Weinberg (1996) model to
show that dwarfs forming well before the peak of reionization ($z>10$) will be
less affected, thus reducing the threshold mass for these objects. This was
modelled in detail by Benitez-Llambay & Frenk (2020), who calculated the
redshift–dependent critical halo mass for star formation before and after
reionization. Assuming instantaneous reionization, they demonstrate a lower
threshold mass before reionization, corresponding to the atomic hydrogen limit
set by a virial temperature of $T_{200}\sim 7000~{}\rm{K}$. After
reionization, the threshold mass is set by the temperature of ionized
hydrogen: $2\times 10^{4}~{}\rm{K}$. Once a halo is in thermal equilibrium
with the background, the metallicity of the gas will not have much effect,
since metal cooling lines are sub-dominant to hydrogen at this temperature.
This critical mass theory connects the total mass growth history to the star
formation history, since halos are only expected to form stars when they are
more massive than the critical mass. If their growth does not keep pace with
the critical mass growth, then they are predicted to cease star formation.
Thus, we can call this a mass growth threshold.

These various models capture the general trends very well and have been
substantiated by multiple galaxy formation simulations (i.e. Oñorbe et al.,
2015; Pereira Wilson et al., 2022). However, the evolution before reionization
is unconstrained observationally, and theoretical predictions span three order
of magnitude for the pre-reionization cut-off mass, covering the range
$10^{5}-10^{8}~{}\,\mathrm{M}_{\odot}$. The lower end of this range is
populated by models assuming that molecular hydrogen cooling leads to
Population III star formation (Machacek et al., 2001; Abel et al., 2002;
Yoshida et al., 2003; Skinner & Wise, 2020; Kulkarni et al., 2021). The upper
end of the range is a result of assuming that atomic hydrogen is the dominant
cooling channel, which happens when the virial temperature rises above
$T_{200}=7000~{}\mathrm{K}$ (e.g., Tegmark et al., 1997; Benitez-Llambay &
Frenk, 2020). This translates to a mass threshold of around
$10^{8}\,\,\mathrm{M}_{\odot}$ at $z=8$.

However, since the critical mass is expected to be redshift–dependent and have
a discontinuity at reionization, the assumed timing of reionization becomes a
decisive parameter here as well. For example, the Tegmark et al. (1997) model
allows for SF in halos of $M_{\mathrm{halo}}=2\times
10^{6}~{}\,\mathrm{M}_{\odot}$ at $z=30$. But this rises steeply to
$10^{8}\,\mathrm{M}_{\odot}$ by $z=9$. New James Webb Space Telescope (JWST)
detections of galaxies at $z\sim 9-17$ show high rest-frame UV luminosities
(Harikane et al., 2022). The authors conclude that this can either be
explained by a higher than expected star formation rate surface density or a
top-heavy stellar initial mass function (IMF). The former could indicate that
the star formation efficiency is higher pre-reionization because it is not
suppressed by a UV background. The latter could indicate that Population III
stars are indeed more massive on average than later stars.

The aim of most of the models presented above is to predict the first
primordial collapse threshold. This means that metal line and metal molecular
cooling are not taken into account. So while this is an extremely interesting
and important question, it really only applies for the very first, isolated
and metal–free halos. However, the average $10^{8}\,\,\mathrm{M}_{\odot}$ halo
at $z=10$ is already expected to be the product of the assembly (merging) of
many smaller halos. Thus, the SFH we would like to understand is actually a
cumulative history of the SF in each of these individual smaller halos. To
complicate matters, the evolution and SF of such a halo at high redshift is
quickly influenced by metal production (i.e. Curti et al., 2022).

After the first (metal–free) Population III (PopIII) stars form, it is
expected that their metal production will enrich the surrounding gas and
contribute to further cooling (i.e. Wise et al., 2012). As shown in Jeon et
al. (2017) and Gutcke et al. (2022), SNe from star forming halos can eject
metals beyond their own virial radii and enrich small neighboring halos that
otherwise would not be able to cool efficiently via either molecular hydrogen
or hydrogen line cooling. Once enriched, the cooling in these small systems is
boosted by metal lines. This process allows certain halos in the vicinity of
larger neighbors with masses down to $\lesssim 10^{3}~{}\,\mathrm{M}_{\odot}$
to host stars prior to reionization. Thus, the assembly of a single
$>10^{8}\,\,\mathrm{M}_{\odot}$ halo can include a combination of primordial
and pre-enriched halos.

Figure 1: Mock V-band surface brightness maps for each simulation at $z=0$.
Each image is $1.4~{}\mathrm{kpc}\times 1.4$ kpc. Notice the colorbar is
rescaled for each image. The simulation is not rotated, only re-positioned to
place the halo in the center.

Lastly, as the halo mass increases and the UVB decreases, some or most of the
metals can be retained within the halo. Combined, this creates favorable
conditions for gas accretion, gas cooling and condensation. When the gas
finally reaches the critical density for self-shielding, SF quickly reignites.

In this work, we aim to study the SFHs of the smallest dwarf galaxies that
survive reionization and place their behavior into the theoretical context of
the mass threshold models. While these models are generally valid, the
detailed behavior of individual systems seems to be sensitive to various
additional processes. In our analysis, we will attempt to disentangle the
causes for the intermittent SF activity displayed by our simulations.
Situations where multiple simulations demonstrate similar timing point to a
cosmological cause that could be interpreted in terms of a larger-scale
synchronicity.

The paper is organized as follows: in Sec. 2 we summarize the most relevant
aspects of the simulation model. In Sec. 3.1 we present the star formation
histories of the five simulations. Then in Sec 3.2 we discuss how these
histories came about through an analysis of the halos’ early evolution. In
Sec. 3.3, we compare the $z=0$ attributes of the simulations with observations
of Local Group dwarf galaxies. Lastly, in Sec. 4 and 5 we discuss our results
in the context of other simulation work and present our main findings.

Table 1: Summary of simulation attributes: virial mass at $z=0$, stellar mass within $R_{200}$ at $z=0$, virial radius, stellar half mass radius, concentration parameter, bulk velocity of the halo, mean SFR until $z=4$, fraction of stars formed in the first 1 Gyr, V-band magnitude at $z=0$, luminosity at $z=0$, and mean stellar metallicity at $z=0$. Name | $M_{200}$ | $M_{\star}$ | $R_{200}$ | $R_{1/2}$ | $c_{\mathrm{max}}$ | $v_{\mathrm{bulk}}$ | $\log\langle\dot{M}\rangle$ | $f_{\mathrm{1Gyr}}$ | $M_{V}$ | $\log L$ | $\log Z$
---|---|---|---|---|---|---|---|---|---|---|---
| $[10^{9}\mathrm{M}_{\odot}]$ | $[10^{6}\mathrm{M}_{\odot}]$ | $[\mathrm{kpc}]$ | $[\mathrm{pc}]$ |  | $[\mathrm{km~{}s}^{-1}]$ | $[\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}]$ |  |  | $[\mathrm{L}_{\odot}]$ | $[\mathrm{Z}_{\odot}]$
HaloA | 3.34 | 9.42 | 28.7 | 311.5 | 13.5 | 177.5 | -2.3 | 0.519 | -11.5 | 6.69 | -1.02
HaloB | 2.82 | 1.26 | 29.8 | 483.2 | 12.7 | 146.6 | -3.2 | 0.992 | -8.8 | 5.79 | -2.60
HaloC | 2.16 | 1.05 | 27.2 | 878.4 | 13.0 | 46.3 | -3.1 | 1.000 | -8.6 | 5.57 | -2.22
HaloD | 1.72 | 1.59 | 25.3 | 119.6 | 24.8 | 120.0 | -2.8 | 0.830 | -9.0 | 5.65 | -1.47
HaloE | 0.29 | 0.22 | 7.2 | 107.0 | 14.7 | 11.0 | -3.7 | 1.000 | -8.4 | 5.55 | -2.52

##  2 Model & Method

In the current study, we present the first sample of five cosmological
simulations run to $z=0$ using the LYRA galaxy formation model (Gutcke et al.,
2021, 2022). LYRA is a comprehensive numerical model that includes a resolved
ISM with a cooling prescription down to $10$K, individual (star-by-star) star
formation, resolved supernovae and a subgrid model for PopIII star enrichment
at high redshift. These prescriptions are implemented in the cosmological,
hydrodynamical moving-mesh code AREPO (Springel, 2010; Pakmor et al., 2016;
Weinberger et al., 2019). We direct the reader to these papers for a detailed
description of the model and the code characteristcs. Below, we will summarize
the most relevant aspects for this study.

The ultra-violet background (UVB) is implemented as an isotropic heating term
following the rates presented in Faucher-Giguère (2020). HI reionization is
set to $z=7.8$, but is $>90\%$ complete by $z=7.3$. HeII reionization is
$>90\%$ complete at $z=3.0$. Gas is self-shielded from this heating following
the fitting function presented in Rahmati et al. (2013). In practice, the
self-shielding factor depends on the gas density and the redshift. The
majority of gas will be fully self-shielded above ${n}_{\rm H}\gtrsim
10^{-2}~{}{\rm cm}^{-3}$. To enable the compatibility of heating and cooling
rates, we choose to use the metal and low temperature cooling rates tabulated
in Ploeckinger & Schaye (2020, “UVB dust1 CR0 G0 shield1”). These were
computed using Cloudy (Ferland et al., 2017) and use the Faucher-Giguère
(2020) UVB as an input. Importantly, instead of the common “1-zone” model,
these tables account for the full absorption of the incident radiation down to
a fixed column density. We include these cooling rates from $z=20$ at the
lowest value tabulated. The simulations are started with zero metallicity,
thus no metal cooling can occur initially. Using the “$10^{6}$” case presented
in Gutcke et al. (2022) as our fiducial model, we implement the enrichment by
PopIII stars in a sub-grid manner. Running the halo finder Subfind (Springel
et al., 2001) on the fly, the metallicity of the gas within the virial radius
is increased from $0$ to $10^{-4}\,\mathrm{Z}_{\odot}$ when a halo crosses the
mass threshold of $M_{\mathrm{PopIII}}=10^{6}\,\mathrm{M}_{\odot}$. Once the
gas metallicity is $10^{-4}\,\mathrm{Z}_{\odot}$ or higher, is it allowed to
cool and form stars.

Star formation follows a Schmidt relation (Schmidt, 1959):

$\dot{M}_{\star}=\varepsilon_{\mathrm{SF}}\frac{M_{\mathrm{gas}}}{t_{\mathrm{ff}}},$
(1)

where $\varepsilon_{\mathrm{SF}}$ is the star formation efficiency parameter
and $t_{\mathrm{ff}}=\sqrt{3\pi/(32G\rho)}$ is the free-fall time of the gas
cell with density $\rho$, and $G$ is the gravitational constant. Gas cells
with $T/\mathrm{K}<100$ and $3<\log({n}_{\rm H}/{\rm cm}^{-3})<4$ are eligible
to form stars at an efficiency of $\varepsilon_{\mathrm{SF}}=2\%$. When
$\log({n}_{\rm H}/{\rm cm}^{-3})>4$, the efficiency is increased to
$\varepsilon_{\mathrm{SF}}=100\%$.

The initial conditions were taken from the EAGLE simulation (Schaye et al.,
2015) following the method presented in Jenkins (2013). The final zoom-in
initial conditions simulate the entire EAGLE box ($L=100\,h^{-1}$ Mpc) at low
resolution. Within this box, we define a high resolution region which
encompasses the complete Lagrangian region of a single dwarf galaxy. Within
this region gas can be refined and de-refined, cool radiatively and form
stars. The DM mass resolution is $\sim 80\,\,\mathrm{M}_{\odot}$, while the
target gas cell mass is set to $4\,\,\mathrm{M}_{\odot}$. The gravitational
softening length is $10$ pc for DM and $4$ pc for the gas and stars. The
galaxies were chosen to fulfill an isolation criterion, namely to not interact
with a larger galaxy for the duration of their lifetime.

##  3 Results

###  3.1 Star formation histories

In this work, we present a sample of five halos in the mass range $1-4\times
10^{9}~{}\,\mathrm{M}_{\odot}$. Fig. 1 shows stellar surface brightness maps
at $z=0$ for each of the five simulations. The general morphology is that of
dwarf spheroidals (dSph), with a smooth extended distribution. However, there
are two additional features of note in these mock images. Firstly, the center
of the main spheroid seems to show signs of nucleation, a central stellar
overdensity distinct from the extended profile. Nucleation is seen in dwarf
galaxies in both the Fornax cluster (Muñoz et al., 2015; Eigenthaler et al.,
2018, Next Generation Fornax Survey) and the Virgo cluster (Grant et al.,
2005, Isaac Newton Telescope Virgo Survey). Secondly, there is a clear
indication of luminous substructure associated with each halo. We summarize
the general attributes of the five simulations in Tab. 1.

In Fig. 2 we present the five star formation histories as a function of
redshift. More precisely, this is an initial–mass–weighted histogram of
formation times of all stars within $R_{200}$ at $z=0$. It includes accreted
stars, but does not include stellar mass loss. The thin lines were computed
using a bin width of $10$ Myr ($500$ Myr for the thick lines). The halos are
ordered from top to bottom by halo mass at $z=0$. HaloA and HaloD show a
similar SFH: they initially quench at $z=7.8$, recommence SF at $z\approx 5$
and continue to form a significant fraction of their stars after reionization.
Nevertheless, it should be pointed out that they are severely affected by
reionization. Firstly, because they initially are quenched for approximately
$400$-$500$ Myr before recommencing SF. Secondly, because the level of SF
after reionization steadily decreases to values around
$10^{-6}-10^{-4}~{}\,\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$, never again
attaining the early (pre-reionization) value of $\sim
10^{-2}~{}\,\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$.

HaloB is also initially quenched by reionization, but has a single burst
around $z=5$, which reduces the central density such that no more stars are
able form at that time. However, 6 billion years later the galaxy rejuvenates
at around $z=0.3$ at a rate of $\sim
10^{-5}~{}\,\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$. HaloC and HaloE are
consistent with being reionization relics. Their star formation ceases at
$z=7.8$, after which they are never again able to assemble sufficient gas to
recommence forming stars. Interestingly, these two galaxies also form their
first stars later than the other three galaxies. While HaloA, HaloB and HaloD
begin forming stars shortly after $z=20$, these two halos do not host stars
until around $z=16$. We also note that the $z=0$ halo mass is neither a good
predictor of the final stellar mass, nor of the shape of the SFH.

###  3.2 Why do certain halos survive Reionization?

In the following, we will investigate the cause of the quenching at
reionization in more detail. Most importantly, we will focus on the reasons
why two (three, if we include HaloB) of the halos recommence their star
forming activity, while the others do not.

Name | $z_{\mathrm{rejuv}}$ | $t_{\mathrm{rejuv}}~{}[\mathrm{Gyr}]$ | $\Delta t~{}[\mathrm{Myr}]$ | $t_{\mathrm{sc}}~{}[\mathrm{Myr}]$
---|---|---|---|---
HaloA | 5.14 | 1.140 | 470 | 155
HaloD | 5.04 | 1.168 | 493 | 163
HaloB | 4.59 | 1.312 | 640 | 128
Table 2: Initial quenching at $z_{\mathrm{qu}}\sim 7.7$ or
$t_{\mathrm{qu}}\sim 0.67~{}\mathrm{Gyr}$ for all halos. Columns are: redshift
of first rejuvenation after quenching, time of rejuvenation, duration of
quenched episode, and mean sound crossing time of reionization–heated halo gas
during the quenched episode.

Figure 2: SFH in bins of 10 Myr (thin lines) and 500 Myr (thick lines). On the
left the $x$-axis is redshift, while on the right it is lookback time in Gyr.
HaloA and HaloD restart SF around $z=5$. HaloB and HaloE remain quenched after
reionization. The grey bands in each panel indicate the time of hydrogen
reionization.

Figure 3: Each of the four figures shows three panels as a function of time
for one simulation, respectively: the top panels show gas density within 40 pc
(solid colored) and within $R_{200}$ (black line). The central panels show the
cooling time measured within the virial radius (colored), the dynamical time
of the halo (black), and the Hubble time (dotted). The lower panels show the
net mass flow rate (outflow rate minus inflow rate) across the spherical shell
with radius $R_{200}$ and width $0.1~{}R_{200}$. The grey bands in each panel
indicate when SF is ongoing. The black vertical dashed line shows the time of
quenching ($z=7.7$), which is just after the central reionization redshift
($z=7.8$).

Fig. 3 presents the evolutionary behavior of the four more massive
simulations, HaloA-D, in the time interval $t=0.2-1.7~{}\mathrm{Gyr}$. Each
subplot shows data from one simulation. The top panel is the average hydrogen
number density within the central $40$ pc (colored line) and within $R_{200}$
(black line). The middle panel is the cooling time of the entire halo (colored
line) computed as

$t_{\mathrm{cool}}=\frac{u}{{\rm d}u/{\rm d}t},$ (2)

where $u$ is the total internal energy of gas cells within the halo. We also
show the dynamical time of the halo (black line), which is

$t_{\mathrm{dyn}}=\left(\frac{3\pi}{32\,G\rho}\right)^{1/2},$ (3)

where $\rho$ is the average total mass density within the halo. Finally, the
bottom panel is the net gas flow rate across a spherical shell centered on
$R_{200}$ with a width of 0.1 $R_{200}$. The grey bands indicate ongoing SF
within the halo at that time. We note that the grey bands differ slightly from
the SFHs in Fig. 2, since the SFHs include all stars that reside within the
main halo at $z=0$. Due to late accretion of luminous substructure, some of
these stars are not yet present in the halo at higher redshifts.

Before reionization, all four simulations show similar behavior: the cooling
time is well below the dynamical time and SF is fairly continuous. The central
gas density varies but has peaks above $1\,{\rm cm}^{-3}$. The outflow and
inflow rates are generally well matched, meaning that SF is self-regulated:
external accretion fuels SF and SNe drive outflows without preventing
accretion. Thus, there is neither run–away cooling nor quenching. This self-
regulation is strongly affected by the onset of reionization.

Hydrogen reionization begins at $z=8.3$ and completes at $z=7.3$, thus
occurring at a mean redshift of $z=7.8$ (black dashed line). During this time,
all gas that is not self-shielded is heated to approximately $1.5\times
10^{4}~{}\mathrm{K}$. In response to this sudden heating, unshielded gas in
the ISM and the entire circum-galactic medium (CGM) must expand. As Shapiro et
al. (2004) show, even the self–shielded gas can eventually become
photo–evaporated. The expansion causes outflows, so there is less net
accretion onto the central region. SF comes to a sudden halt in all halos due
to a lack of fuel. Cooling times of the gas increase above the dynamical time
and even above the Hubble time (grey dotted line). The net flow rate may be
outward due to the last SNe driving outflows into the now ionized medium,
making the central gas density drop drastically.

Let us assume the majority of gas is in hydrostatic equilibrium before the
onset of reionization. The gas is neither isothermal nor in thermal
equilibrium with an external radiation field (since there is none in this
model). A pressure (and temperature) gradient across the halo supports the
non-SF gas. We can approximate reionization as an instantaneous (non-
adiabatic) heating process.

The sudden heating of the gas by reionization raises the pressure. To
reestablish equilibrium, the gas in the halo will expand on the sound crossing
time scale, which is the time it takes to cross the halo when travelling at
the local sound speed. It is straightforward to compute the adiabatic sound
speed assuming an ideal gas:

$c_{s}=\sqrt{\gamma\frac{k_{\mathrm{B}}T}{\mu m_{p}}}\approx
18.4~{}\mathrm{km~{}s}^{-1}\left(\frac{T}{1.5\times
10^{4}~{}\mathrm{K}}\right)^{1/2},$ (4)

where $\gamma=5/3$ is the adiabatic index, $k_{\mathrm{B}}$ is the Boltzmann
factor, $T$ is the temperature, $m_{p}$ is the proton mass, and we assume a
mean molecular weight of ionized hydrogen and singly ionized helium after
reionization, $\mu=2/(3X+1)\approx 0.61$, where $X=0.76$ is the primordial
hydrogen mass fraction. The two halos with early rejuvenation, HaloA and HaloD
undergo the same quenching and recommencement of SF at approximately the same
time. They have an average virial radius of $R_{200}\approx 3~{}\mathrm{kpc}$
between $z=7.8-5$. Thus, the sound crossing time of these halos is
$t_{\mathrm{sc}}\approx 160$ Myr.

Due to the sudden heating and subsequent expansion, some gas may become
unbound, as seen in the net outflow rates. However, most of the gas will cool
and re-collapse due to gravity. The net flow rates reverse and inflow begins
to dominate after around one sound crossing time, driving gas back inside the
virial radius. The initial expansions and subsequent recollapse can be seen as
oscillation (or “breathing”) of the halo in response to reionization.

The total gas density (black line) in HaloB and HaloC steadily decreases
immediately after reionization. There is photo-evaporation of the gas at and
beyond the virial radius. Accretion is halted and outflows dominate. However,
within one sound crossing time all four halos reestablish a central gas
reservoir. The minimum density required to do this is the self-shielding
threshold of $\sim 10^{-2}~{}{\rm cm}^{-3}$. After that, the central gas
densities (colored lines) either remain constant or rise steadily. This
accretion occurs despite the total halo density (black lines) remaining
constant or decreasing. This implies that most gas does not escape the halo,
but is initially suspended in the CGM until it is slowly re-accreted onto the
central zone.

As the central densities increase, the cooling time decreases. Once it drops
below the dynamical time again, SF can start again. This happens after
$\sim$500 Myr in HaloA, HaloB and HaloD. The time required for this is a
combination of the sound crossing time and the time required for outwards
angular momentum transport and radiative cooling. We suggest that this
$\sim$500 Myr break results from intrinsic properties of these halos that are
at the edge between quenching and survival. The fact that all three halos
display the same behavior can be interpreted as synchronicity with a
cosmological (reionization) cause.

Figure 4: Top: Evolution with redshift of the virial temperature of the halos.
The horizontal dashed grey line shows $T=10^{4}\,\mathrm{K}$, the limit for
atomic cooling in halos. The vertical dashed grey line shows the approximate
timing of reionization when all halos halt their star formation. For
comparison, we also plot the IGM temperature. Bottom: Evolution with redshift
of the virial mass of the halos. The opacity of the lines indicate when SF is
ongoing (full) or quenched (opaque). For comparison, we plot the model of
Benitez-Llambay & Frenk (2020) as a solid black line. The horizontal grey
dashed line indicates the PopIII enrichment threshold in our model (before
reionization) while we show the atomic cooling threshold of $10^{4}$ K for
$z<7.7$. The vertical grey dashed line denotes the quenching redshift of
$z=7.7$.

HaloB shows a slightly different behavior than HaloA and HaloD after
reionization. Already before reionization it shows a more bursty and sporadic
SF with periods of lower gas densities. There is a final burst of SF at the
time of reionization and some subsequent net outflow. However, the net flow
does not turn around to become an inflow in the same way as in the other two
halos. Instead, there is a low but persistent net outflow, meaning that the
halo is photo-evaporating at its boundary. So while we can see a slow increase
in the central gas density, this is material from the CGM, not new accretion.
The central density reaches a sufficiently high value that SF recommences but
is interrupted after a single burst, sending the cooling time to values
greater than the Hubble time. SF is quenched for the next 6 Gyr. Then, at
$z\approx 0.3$, it recommences at a low rate.

Lastly, we will look at HaloC (which is very similar to HaloE and not shown
here for simplicity). SF is quenched slightly before the peak of reionization,
showing that the halo is especially sensitive to the external radiation field.
The cooling time increases rapidly and the net flow rates are very small (note
the $y$-axis is rescaled), showing no large-scale accretion. Interestingly,
while this halo does not manage to form stars again, there is nevertheless an
indication of activity around the same time that the other halos rejuvenate.
The cooling time decreases below the Hubble time and the gas density rises
slightly. After this, the activity ceases.

Figure 5: Mass loading factor as a function of time, using the mean SFR until
$z=4$ as opposed to the instantaneous SFR.

HaloA-D were chosen to have $z=0$ halo masses within a factor of two of each
other. Nevertheless, their growth histories vary. In Fig. 4 we show the virial
temperature of each halo as a function of redshift:

$T_{200}=\frac{\mu m_{p}}{2k_{\mathrm{B}}}\frac{GM_{200}}{R_{200}},$ (5)

where $M_{200}$ is the viral mass and $R_{200}$ is the virial radius. To guide
the eye, we have also plotted the atomic cooling threshold of $10^{4}$ K and
the quenching redshift of $z=7.7$. The magenta line shows the IGM temperature,
computed following Puchwein et al. (2015), namely the median gas temperature
of IGM gas within 5% of the mean cosmic density. The IGM temperature begins to
rise at $z=8.3$ according to our UVB, and reaches $10^{4}~{}\mathrm{K}$ by
$z=7.3$. Starting at $z=4$ we begin to see HeII reionization heating the IGM.
After $z=2$ the IGM temperature drops below $10^{4}~{}\mathrm{K}$ due to the
expansion of the universe.

We see that HaloD displays the fastest early growth, rising above $10^{4}$ K
as early as $z\sim 10$. This is also seen in the concentration parameter (see
Tab. 1). The concentration is computed following Bullock et al. (2001)
assuming an NFW profile shape (Navarro et al., 1997) and using the
$r_{\mathrm{max}}$–method:

$c=2.16\frac{R_{200}}{r_{\mathrm{max}}},$ (6)

where $r_{\mathrm{max}}$ is the radius at the peak of the rotation curve,
$v_{\mathrm{max}}$. Large concentration parameters are known to be indicative
of early growth (e.g. Croton et al., 2007), since the central halo builds up
early, when the Universe was denser. The fact that this affects halo
properties, such as SF, is known as assembly bias. HaloA also grows
sufficiently to pass this crucial threshold before the onset of reionization,
despite having a much lower concentration. HaloB and HaloC cross this
threshold after reionization, while HaloE never grows sufficiently.

Fig. 5 shows the mass loading factor as a function of redshift for all five
simulations. The mass loading factor is the ratio of the outflow rate over the
mean SFR, where we have measured the mean value from the time of the first
star to $z=4$:

$\eta_{M}=\frac{\dot{M}_{\mathrm{out}}}{\langle\dot{M}_{\star}\rangle}.$ (7)

Since all simulations have discontinuous SFHs, the mean SFR is lower than if
we had only accounted for times of ongoing SF. Therefore, the mass loading
factors shown are lower than the instantaneous values would be. The
time–averaged mass loading factors from the time the first star forms to $z=4$
for the simulations in alphabetical order are: $33.7$, $26.3$, $3.7$, $6.5$,
and $7.0$. It is interesting to note that the mean values are above unity in
all cases. This means the halos are losing substantial gas mass throughout
their evolution. We show this even more clearly in Fig. 6, in terms of the
baryon fraction across redshift. Until $z\approx 15$, all halos are consistent
with hosting baryons at the cosmic baryon fraction (grey dashed line). Once SF
sets in, mass loss occurs. While certain halos are able to regain higher
fractions for some time, the final values tend to assume extremely low values
of $1-10\%$ of the mean cosmic value. It is of interest to note that HaloE
does not form new stars after $z\approx 7.7$ but nevertheless shows a mass
loading above unity and a decreasing baryon fraction. This must therefore be
attributed to photo-evaporation due to the UVB.

Figure 6: Baryon fraction normalized to the cosmic baryon fraction,
$\Omega_{\mathrm{b}}/\Omega_{\mathrm{m}}$, as a function of redshift for each
simulation. The vertical grey band indicates the time of hydrogen
reionization.

###  3.3 What are the z=0 properties?

Figure 7: The peak of the stellar rotation curve, $v_{\mathrm{rot}}$, the
1D-velocity dispersion within the half-mass radius, $\sigma_{\star}$ and their
ratio as a function of total stellar mass. Grey data points with error bars
show the observations from Wheeler et al. (2017). HaloA and HaloD are
rotationally supported.

Now that we have discussed the SFHs and the reason why certain halos survive
reionization, let us turn to the $z=0$ properties. Do these properties conform
to Local Group observations? Do they betray their different histories?

Firstly, in Fig. 7 we look at the kinematics of the stars at $z=0$. In the
left panel, we show the peak value of the stellar rotation curve. The central
panel is the stellar velocity dispersion measured within the half–mass radius.
To be precise, this is the average of the three independent ($x$, $y$ and $z$)
line-of-sight velocity dispersions. The right panel shows the ratio of these
two values, giving us an indication of the dominant mode of support against
gravitational collapse. In each panel, the grey data points are observations
presented in Wheeler et al. (2017). The rotation velocity of our five halos as
a function of total stellar mass matches the data well. There is a distinct
increase in velocity with increasing stellar mass which reflects the same
trend seen in the data. Also, the velocity dispersions in the central panel
match the data. However, they seem to be slightly low in the simulations. This
is likely caused by the “softened” treatment of gravitational N–body forces
below the softening length, $\varepsilon_{\star}$, which we have set to
$4\,{\rm pc}$ in this work. This means that the gravitational acceleration
acting on stars passing within $\varepsilon_{\star}$ of another particle is
reduced, a choice made for computational efficiency. This can affect the
estimated value of the velocity dispersion when interactions such as
encounters closer than the softening length and binary stars interacting with
a third star occur. Such interactions would cause an exchange of velocities,
and, therefore, increase the total dispersion of the system. Since they are
not accounted for here, the velocity dispersion is on the low end of what is
expected from observations. Finally, in the right panel we see that the four
lower mass systems are dispersion-dominated (values below unity). HaloA, on
the other hand, sits close to unity, consistent with the appearance of more
rotation-dominated systems at these stellar masses in the observational data.

Figure 8: Metallicity as a function of V-band magnitude. Grey data points show
LG dwarf observations from McConnachie (2012). Figure 9: Age–metallicity
relation for HaloA. The colorbar shows the logarithm of the mean SFR in each
bin of age and metallicity. Age bins are 500 Myr wide, while metallicity bins
are 0.25 dex. Compare this figure to Fig. 8a of Savino et al. (2019) showing
the age-metallicity relation of Tucana dSph.

Fig. 8 shows the magnitude–metallicity relation for our simulations at $z=0$
and the LG data from McConnachie (2012). The metallicity is the log–averaged
$z=0$ metallicity of all stars within $R_{200}$. Our simulations span the
range of metallicities in the data, between $\log Z/\,\mathrm{Z}_{\odot}=-2.5$
to $-1$. As Agertz et al. (2020) showed, the metallicity relation is a useful
indicator for the strength of feedback. Metallicities above the relation
indicate a too weak feedback prescription, since metals are not ejected out of
the ISM, instead accumulating and forming high metallicity stars. On the other
hand, if the feedback is too strong then metals are preferentially ejected
with the SN bursts and impoverish the ISM. The fact that this work spans the
data is an encouraging indication that resolving the cooling radius for each
SN self-consistently calibrates the strength of feedback.

It is unsurprising that the two halos that survive reionization also show the
highest metallicities, since later SF produces more enriched stars. To show
this in detail, we plot the age–metallicity diagram for HaloA in Fig. 9. SF
starts before reionization at the lowest metallicity allowed for SF, $\log
Z/\,\mathrm{Z}_{\odot}=-4$. Within the first Gyr, the metallicity increases to
$>-2$. After $z=4$, no new stars form below $\log
Z/\,\mathrm{Z}_{\odot}=-1.5$. However, the mean SFR is also reduced
significantly. Recent Early Release Observations with the James Webb Space
Telescope (JWST) using the Near-Infrared Spectrograph (NIRSpec) indicate a
rapid build-up of metallicity at $z>7$ (Arellano-Córdova et al., 2022). The
metallicity build-up within our halos prior to the onset of reionization is
also rapid, with the first solar metallicity stars forming as early as $z>10$.

Figure 10: V-band magnitude as a function of stellar half mass (light) radius,
$R_{1/2}$. Grey data points show LG dwarf observations from McConnachie
(2012). Filled diamonds shows half mass radii. Small back diamonds also show
the Milky Way globular clusters from Harris (1996). The grey shaded area are
sizes below the dark matter softening length of our simulations. The diagonal
light blue lines are lines of constant surface brightness. Figure 11:
Cumulative radial distribution of stars within $R=1~{}\mathrm{kpc}$ for HaloA.
The red line shows the distribution for stars formed before reionization,
while the blue line is all stars formed afterwards. Clearly, the younger stars
are more centrally concentrated.

Finally, we also show the size–magnitude relation in Fig. 10. Colored diamonds
give the stellar half–mass radii in the simulations. Again, the grey data
points show LG dwarfs from McConnachie (2012). For comparison, black diamonds
show Milky Way globular clusters (GCs) from Harris (1996, 2010 edition).
Clearly, the simulated sizes are within expectations of the data. All five
simulations are well placed within the classical dwarf regime, not in the
ultra-faint regime ($M_{V}\lesssim-6$). The two halos with more recent SF have
smaller sizes. This is due to the fact that younger stars, which are brighter
and have not lost mass, are more centrally located. We show this in Fig. 11
for HaloA. This is consistent with Fig. 11b of Savino et al. (2019) for Tucana
dSph. In Fig. 12 we see that the central 0.5 kpc of HaloA are also dominated
by a rotating disk, which is composed primarily of younger stars. This is line
with 7, in which only HaloA shows signs of rotational support.

Figure 12: Line-of-sight ($z$-axis) stellar velocity (left), velocity
dispersion (center) and skewness (right) of the central $5R^{\star}_{1/2}$ for
HaloA. The halo was rotated such that the angular momentum vector of the stars
points in the $y$-direction, thus we see a rotating disk edge-on. Some stellar
substructure is visible.

##  4 Discussion

Cosmological simulations at mass resolutions of $<10~{}\,\mathrm{M}_{\odot}$
are not yet very common. But this level of resolution is crucial for the
galaxy mass scale we study here, since it allows us to resolve the multi-phase
ISM and individual SN bursts without parameters tuned to match observations.
We can compare our results with the EDGE simulations (Agertz et al., 2020)
that have a mass resolution of $20\,\,\mathrm{M}_{\odot}$. In Rey et al.
(2020, hereafter R20), the SFHs of four galaxies from this project were
presented in detail. The authors show that reionization initially quenches all
four. However, two are able to reignite SF at low rates after $z=1$. These two
cases are reminiscent of HaloB in our sample. Of the two quenched halos in
R20, the “gas-rich” simulation is similar to HaloC presented here. While
neither forms stars after reionization, they are very close to doing so.

R20 show that a small “genetic” modification that increases the final mass can
incite this halo to rejuvenate. In our work, HaloC is the same initial
condition as presented in Gutcke et al. (2022). In that paper, we showed that
small modifications to the treatment of metals at high redshift (for example
setting $M_{\mathrm{PopIII}}=10^{7}\,\,\mathrm{M}_{\odot}$) also had the
effect of inducing rejuvenation. Thus, both of these halos seem to demonstrate
properties at the extreme edge between survival and quenching. Minor model
variations are able to tip the evolution in one direction or another. Both
halos have a very similar $z=0$ total masses of $2.5\times
10^{9}~{}\,\mathrm{M}_{\odot}$ (R20) and $2.16\times
10^{9}~{}\,\mathrm{M}_{\odot}$ (HaloC), respectively. One additional point to
note, however, is that the stellar mass of the two are almost a factor of two
different: $5.9\times 10^{5}~{}\,\mathrm{M}_{\odot}$ (R20) and $1.05\times
10^{6}\,\mathrm{M}_{\odot}$ (HaloC).

Wheeler et al. (2019) present FIRE dwarfs with a mass resolution of
$30\,\,\mathrm{M}_{\odot}$. They show consistency with Oñorbe et al. (2015) in
which star formation in ultra-faint dwarf galaxies is quenched by
reionization, but residual self-shielded gas allows stars to continue forming
until $z\sim 2.5$. The lack of fresh accretion, which is caused by
reionization, then prevents SF after $z=2$. They also predict that all
simulations producing classical dwarfs continue to form stars until at least
$z=0.5$.

These results are somewhat in conflict with our work here. This is primarily
caused by quite different stellar masses in the same halo mass ranges. For
example, the “$m09_{30}$” simulation in Wheeler et al. (2019) has a halo mass
of $M_{\mathrm{halo}}=2.5\times 10^{9}~{}\,\mathrm{M}_{\odot}$, but only
produces $1.2\times 10^{4}~{}\,\mathrm{M}_{\odot}$ in stars. For the same halo
mass, this is two orders of magnitude lower than the stellar mass formed in
our work here. This discrepancy is already set before the onset of
reionzation, since the SFR in our work is $\gtrsim
10^{-2}~{}\,\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$, and all five galaxies form
$\gtrsim 10^{5}~{}\,\mathrm{M}_{\odot}$ before $z=10$. “$m09_{30}$” on the
other hand forms stars at a rate of $\lesssim
10^{-4}~{}\,\mathrm{M}_{\odot}~{}\mathrm{yr}^{-1}$ before reionization,
producing an order of magnitude less stars by the same time.

It is plausible that these large differences so early on in the evolution are
caused by a combination of two model choices. For one, the FIRE model does not
account for PopIII enrichment. Possibly more importantly, the stellar feedback
prescription is known to be very strong and bursty in the FIRE model. It is
likely that this suppresses early SF too much in their model. This explanation
is consistent with the extremely low metallicities presented in the same work.
Encouragingly, the updated FIRE-3 model (Hopkins et al., 2022) appears to
present altered predictions for precisely this halo mass range ($1-5\times
10^{9}~{}\,\mathrm{M}_{\odot}$). Fig. 9 of that paper shows stellar masses one
order of magnitude higher in the new model, bringing their results in much
closer agreement to ours. The remaining differences can be attributed to the
additional feedback prescriptions included in FIRE-3, such as photo-ionization
in HII regions, photo-electric heating, and cosmic ray heating. Future work on
the LYRA model will include these and then a more conclusive comparison can be
made.

##  5 Conclusions

We have carried out five ultra high resolution cosmological zoom-in
simulations of dwarf galaxies in the mass range $1-4\times
10^{9}~{}\,\mathrm{M}_{\odot}$. At
$M_{\mathrm{target}}=4~{}\,\mathrm{M}_{\odot}$, these simulations constitute
the highest resolution zoom-in simulations run fully cosmologically to $z=0$
to date. Additionally, they include a resolved multi-phase ISM and
individually resolved SNe explosions. As our galaxies were selected to form in
isolation and distanced from larger neighbors, our simulations probe universal
aspects of cosmological evolution and of the effects of reionization, without
being influenced by additional effects caused by environment or infall.

The halo mass range chosen for our study covers the edge of where galaxy
formation is expected. As such, it is unsurprising that only two out of the
five simulations form a significant fraction of stars after reionization. To
understand why these two survive while the other three remain quenched, we
analyze the star formation histories and the early evolution to $z=4$ in
detail. Our main results are as follows:

  1. i)

In the studied mass range, the $z=0$ halo mass is not a good predictor for the
survival after reionization.

  2. ii)

Framing the evolution and SFH of dwarfs in terms of an evolving mass threshold
is useful in a general sense, but fails to capture the details of individual
galaxies.

  3. iii)

In our sample, halos that grow to $T_{200}>10^{4}~{}\mathrm{K}$ before
reionization form stars after reionization.

  4. iv)

The two surviving halos are initially quenched by reionization for around 500
Myr, before recommencing SF.

  5. v)

The non-adiabatic heating caused by reionization elicits an oscillation or
breathing of the entire gas within the halo. Initially, the gas expands, then
slows and, if still bound, re-collapses onto the halo. This final step can be
followed by a SF episode.

  6. vi)

Late-time ($z\sim 0.3$) reignition of SF occurs in one of our halos. The halo
growth and gas accretion leading up to the reignition are concurrent with the
drop of the strength of the UVB after $z=2$.

To place our simulations into the context of Local Group dwarf galaxy
observations, we also compare them to various $z=0$ properties of our
simulations. Our main findings are:

  1. i)

The feedback strength resulting from resolved SNe produces a
magnitude–metallicity relation as expected from observations.

  2. ii)

Younger stars, when they exist, are more centrally concentrated and in a
rotating disk.

  3. iii)

Younger stars are more metal-rich, leading to a higher mean metallicity.

  4. iv)

The LYRA model predicts that the examined halo mass range produces classical
dwarfs.

Additionally, we would like to stress that this work is based on five
simulations only. So, while this does constitute a small sample, certain
evolutionary behavior may still be purely stochastic. Also, while most of our
model parameters are either set by first principles or are directly motivated
by observations instead of being calibrated for, there remains some room for
modified parameter choices or assumptions. Some of the results presented here
may be sensitive to the details of these choices.

Finally, our goals for the future are to expand the sample of simulations and
the parameter space covered by this work, as well as to include additional
physics, in particular magneto-hydrodynamics and cosmic rays. Thus we hope to
present in forthcoming work refined predictions for dwarf galaxy properties as
a function of halo mass and be able to either confirm or extend common
extrapolations down into the low mass galaxy regime.

## Acknowledgments

We would like to thank Romain Teyssier and Eve C. Ostriker for useful
discussions. TAG acknowledges support by NASA through the NASA Hubble
Fellowship grant $\\#$HF2-51480 awarded by the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., for NASA, under contract NAS5-26555. CP acknowledges
support by the European Research Council under ERC-AdG grant
PICOGAL-101019746. GLB acknowledges support from the NSF (AST-2108470, XSEDE
grant MCA06N030), a NASA TCAN award 80NSSC21K1053, and the Simons Foundation.
TN acknowledges support from the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy – EXC-2094 –
390783311 from the DFG Cluster of Excellence ‘ORIGINS’. We acknowledge the
computing time provided by the Leibniz Rechenzentrum (LRZ) of the Bayrische
Akademie der Wissenschaften on the machine SuperMUC-NG (pn73we). This research
was also carried out on the High Performance Computing resources of the FREYA
and COBRA clusters at the Max Planck Computing and Data Facility (MPCDF,
https://www.mpcdf.mpg.de) in Garching operated by the Max Planck Society
(MPG).

## References

  * Abel et al. (2002) Abel, T., Bryan, G. L., & Norman, M. L. 2002, Science, 295, 93, doi: 10.1126/science.295.5552.93
  * Agertz et al. (2020) Agertz, O., Pontzen, A., Read, J. I., et al. 2020, MNRAS, 491, 1656, doi: 10.1093/mnras/stz3053
  * Arellano-Córdova et al. (2022) Arellano-Córdova, K. Z., Berg, D. A., Chisholm, J., et al. 2022, arXiv e-prints, arXiv:2208.02562.  https://arxiv.org/abs/2208.02562
  * Babul & Rees (1992) Babul, A., & Rees, M. J. 1992, MNRAS, 255, 346, doi: 10.1093/mnras/255.2.346
  * Benitez-Llambay & Frenk (2020) Benitez-Llambay, A., & Frenk, C. 2020, MNRAS, 498, 4887, doi: 10.1093/mnras/staa2698
  * Benson et al. (2002) Benson, A. J., Lacey, C. G., Baugh, C. M., Cole, S., & Frenk, C. S. 2002, MNRAS, 333, 156, doi: 10.1046/j.1365-8711.2002.05387.x
  * Bullock et al. (2001) Bullock, J. S., Kolatt, T. S., Sigad, Y., et al. 2001, MNRAS, 321, 559, doi: 10.1046/j.1365-8711.2001.04068.x
  * Croton et al. (2007) Croton, D. J., Gao, L., & White, S. D. M. 2007, MNRAS, 374, 1303, doi: 10.1111/j.1365-2966.2006.11230.x
  * Curti et al. (2022) Curti, M., D’Eugenio, F., Carniani, S., et al. 2022, arXiv e-prints, arXiv:2207.12375.  https://arxiv.org/abs/2207.12375
  * Dijkstra et al. (2004) Dijkstra, M., Haiman, Z., Rees, M. J., & Weinberg, D. H. 2004, ApJ, 601, 666, doi: 10.1086/380603
  * Efstathiou (1992) Efstathiou, G. 1992, MNRAS, 256, 43P, doi: 10.1093/mnras/256.1.43P
  * Eigenthaler et al. (2018) Eigenthaler, P., Puzia, T. H., Taylor, M. A., et al. 2018, ApJ, 855, 142, doi: 10.3847/1538-4357/aaab60
  * Faucher-Giguère (2020) Faucher-Giguère, C.-A. 2020, MNRAS, 493, 1614, doi: 10.1093/mnras/staa302
  * Ferland et al. (2017) Ferland, G. J., Chatzikos, M., Guzmán, F., et al. 2017, Rev. Mexicana Astron. Astrofis., 53, 385\.  https://arxiv.org/abs/1705.10877
  * Gnedin (2000) Gnedin, N. Y. 2000, ApJ, 542, 535, doi: 10.1086/317042
  * Grant et al. (2005) Grant, N. I., Kuipers, J. A., & Phillipps, S. 2005, MNRAS, 363, 1019, doi: 10.1111/j.1365-2966.2005.09518.x
  * Gutcke et al. (2021) Gutcke, T. A., Pakmor, R., Naab, T., & Springel, V. 2021, MNRAS, 501, 5597, doi: 10.1093/mnras/staa3875
  * Gutcke et al. (2022) —. 2022, MNRAS, 513, 1372, doi: 10.1093/mnras/stac867
  * Hambrick et al. (2009) Hambrick, D. C., Ostriker, J. P., Naab, T., & Johansson, P. H. 2009, ApJ, 705, 1566, doi: 10.1088/0004-637X/705/2/1566
  * Harikane et al. (2022) Harikane, Y., Ouchi, M., Oguri, M., et al. 2022, arXiv e-prints, arXiv:2208.01612.  https://arxiv.org/abs/2208.01612
  * Harris (1996) Harris, W. E. 1996, AJ, 112, 1487, doi: 10.1086/118116
  * Hoeft et al. (2006) Hoeft, M., Yepes, G., Gottlöber, S., & Springel, V. 2006, MNRAS, 371, 401, doi: 10.1111/j.1365-2966.2006.10678.x
  * Hopkins et al. (2022) Hopkins, P. F., Wetzel, A., Wheeler, C., et al. 2022, arXiv e-prints, arXiv:2203.00040.  https://arxiv.org/abs/2203.00040
  * Jenkins (2013) Jenkins, A. 2013, MNRAS, 434, 2094, doi: 10.1093/mnras/stt1154
  * Jeon et al. (2017) Jeon, M., Besla, G., & Bromm, V. 2017, ApJ, 848, 85, doi: 10.3847/1538-4357/aa8c80
  * Kulkarni et al. (2021) Kulkarni, M., Visbal, E., & Bryan, G. L. 2021, ApJ, 917, 40, doi: 10.3847/1538-4357/ac08a3
  * Macciò et al. (2017) Macciò, A. V., Frings, J., Buck, T., et al. 2017, MNRAS, 472, 2356, doi: 10.1093/mnras/stx2048
  * Machacek et al. (2001) Machacek, M. E., Bryan, G. L., & Abel, T. 2001, ApJ, 548, 509, doi: 10.1086/319014
  * McConnachie (2012) McConnachie, A. W. 2012, AJ, 144, 4, doi: 10.1088/0004-6256/144/1/4
  * McQuinn et al. (2010) McQuinn, K. B. W., Skillman, E. D., Cannon, J. M., et al. 2010, ApJ, 721, 297, doi: 10.1088/0004-637X/721/1/297
  * Muñoz et al. (2015) Muñoz, R. P., Eigenthaler, P., Puzia, T. H., et al. 2015, ApJ, 813, L15, doi: 10.1088/2041-8205/813/1/L15
  * Navarro et al. (1997) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493, doi: 10.1086/304888
  * Oñorbe et al. (2015) Oñorbe, J., Boylan-Kolchin, M., Bullock, J. S., et al. 2015, MNRAS, 454, 2092, doi: 10.1093/mnras/stv2072
  * Okamoto et al. (2008) Okamoto, T., Gao, L., & Theuns, T. 2008, MNRAS, 390, 920, doi: 10.1111/j.1365-2966.2008.13830.x
  * Olsen et al. (2021) Olsen, C., Gawiser, E., Iyer, K., et al. 2021, ApJ, 913, 45, doi: 10.3847/1538-4357/abf3c2
  * Pakmor et al. (2016) Pakmor, R., Springel, V., Bauer, A., et al. 2016, MNRAS, 455, 1134, doi: 10.1093/mnras/stv2380
  * Pereira Wilson et al. (2022) Pereira Wilson, M., Navarro, J., Santos Santos, I., & Benitez Llambay, A. 2022, arXiv e-prints, arXiv:2206.05338.  https://arxiv.org/abs/2206.05338
  * Ploeckinger & Schaye (2020) Ploeckinger, S., & Schaye, J. 2020, MNRAS, 497, 4857, doi: 10.1093/mnras/staa2172
  * Puchwein et al. (2015) Puchwein, E., Bolton, J. S., Haehnelt, M. G., et al. 2015, MNRAS, 450, 4081, doi: 10.1093/mnras/stv773
  * Rahmati et al. (2013) Rahmati, A., Pawlik, A. H., Raičević, M., & Schaye, J. 2013, MNRAS, 430, 2427, doi: 10.1093/mnras/stt066
  * Rees (1986) Rees, M. J. 1986, MNRAS, 218, 25P, doi: 10.1093/mnras/218.1.25P
  * Rey et al. (2020) Rey, M. P., Pontzen, A., Agertz, O., et al. 2020, MNRAS, 497, 1508, doi: 10.1093/mnras/staa1640
  * Savino et al. (2019) Savino, A., Tolstoy, E., Salaris, M., Monelli, M., & de Boer, T. J. L. 2019, A&A, 630, A116, doi: 10.1051/0004-6361/201936077
  * Schaye et al. (2015) Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521, doi: 10.1093/mnras/stu2058
  * Schmidt (1959) Schmidt, M. 1959, ApJ, 129, 243, doi: 10.1086/146614
  * Shapiro et al. (2004) Shapiro, P. R., Iliev, I. T., & Raga, A. C. 2004, MNRAS, 348, 753, doi: 10.1111/j.1365-2966.2004.07364.x
  * Skinner & Wise (2020) Skinner, D., & Wise, J. H. 2020, MNRAS, 492, 4386, doi: 10.1093/mnras/staa139
  * Springel (2010) Springel, V. 2010, MNRAS, 401, 791, doi: 10.1111/j.1365-2966.2009.15715.x
  * Springel et al. (2001) Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726, doi: 10.1046/j.1365-8711.2001.04912.x
  * Tegmark et al. (1997) Tegmark, M., Silk, J., Rees, M. J., et al. 1997, ApJ, 474, 1, doi: 10.1086/303434
  * Thoul & Weinberg (1996) Thoul, A. A., & Weinberg, D. H. 1996, ApJ, 465, 608, doi: 10.1086/177446
  * Weinberger et al. (2019) Weinberger, R., Springel, V., & Pakmor, R. 2019, arXiv e-prints, arXiv:1909.04667.  https://arxiv.org/abs/1909.04667
  * Weisz et al. (2014) Weisz, D. R., Dolphin, A. E., Skillman, E. D., et al. 2014, ApJ, 789, 147, doi: 10.1088/0004-637X/789/2/147
  * Wheeler et al. (2017) Wheeler, C., Pace, A. B., Bullock, J. S., et al. 2017, MNRAS, 465, 2420, doi: 10.1093/mnras/stw2583
  * Wheeler et al. (2019) Wheeler, C., Hopkins, P. F., Pace, A. B., et al. 2019, MNRAS, 490, 4447, doi: 10.1093/mnras/stz2887
  * Wise et al. (2012) Wise, J. H., Turk, M. J., Norman, M. L., & Abel, T. 2012, ApJ, 745, 50, doi: 10.1088/0004-637X/745/1/50
  * Yoshida et al. (2003) Yoshida, N., Abel, T., Hernquist, L., & Sugiyama, N. 2003, ApJ, 592, 645, doi: 10.1086/375810


††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Resonant transport and near-field effects in photonic glasses

Geoffroy J. Aubry  Fachbereich Physik, Universität Konstanz, 78457 Konstanz,
Germany  Lukas Schertel  Fachbereich Physik, Universität Konstanz, 78457
Konstanz, Germany  Physik-Institut, Universität Zürich, 8057 Zürich,
Switzerland  Mengdi Chen  Fachbereich Chemie, Universität Konstanz, 78457
Konstanz, Germany  Henrik Weyer  Fachbereich Physik, Universität Konstanz,
78457 Konstanz, Germany  Christof M. Aegerter  Physik-Institut, Universität
Zürich, 8057 Zürich, Switzerland  Sebastian Polarz  Fachbereich Chemie,
Universität Konstanz, 78457 Konstanz, Germany  Helmut Cölfen  Fachbereich
Chemie, Universität Konstanz, 78457 Konstanz, Germany  Georg Maret
Fachbereich Physik, Universität Konstanz, 78457 Konstanz, Germany

(August 28, 2024)

###### Abstract

A fundamental quantity in multiple scattering is the transport mean free path
the inverse of which describes the scattering strength of a sample. In this
paper, we emphasize the importance of an appropriate description of the
effective refractive index $n_{\mathrm{eff}}$ in multiple light scattering to
accurately describe the light transport in dense photonic glasses. Using
$n_{\mathrm{eff}}$ as calculated by the energy-density coherent potential
approximation we are able to predict the transport mean free path of
monodisperse photonic glasses. This model without any fit parameter is in
qualitative agreement with numerical simulations and in fair quantitative
agreement with spectrally resolved coherent backscattering measurements on new
specially synthesized polystyrene photonic glasses. These materials exhibit
resonant light scattering perturbed by strong near-field coupling, all
captured within the model. Our model might be used to maximize the scattering
strength of high index photonic glasses, which are a key in the search for
Anderson localization of light in three dimensions.

Transport phenomena are omnipresent in nature, governing many processes in
chemistry, biology, physics, and engineering. Systems as diverse as electrons
Abrahams  _et al._ (1979) and ultrasound Bayer and Niederdränk (1993) in
condensed matter, mechanical waves in the earth Larose  _et al._ (2004), cold
atoms in an optical trap Jendrzejewski  _et al._ (2012), and light in
disordered photonic materials Wolf and Maret (1985); *Albada1985 share the
same physical principle Akkermans and Montambaux (2007); Lagendijk  _et al._
(2009). Optical experiments are especially appealing because of the absence of
photon-photon interaction (unlike in electronic systems) and the existence of
relatively high index scattering media such as photonic crystals and glasses.
Moreover, optical transport experiments have reached an unprecedented accuracy
thanks to the great technological development of sources (e.g., lasers),
detectors (e.g., CCDs), and time resolution. All these progresses allow the
realization of table-top experiments which highlight the richness of transport
phenomena.

Wave transport in a diluted disordered suspension of scatterers can be
described by the sole far-field properties of the single scatterers. On
increasing concentration, however, interference effects due to scatterer-
scatterer position correlation need to be taken into account Fraden and Maret
(1990). In this description, the scattering cross section is still the single
scatterer one, which, in general, is calculated in the far field. In optics,
this approach is expected to fail as soon as the photon scattering mean free
path $\ell_{\mathrm{s}}$ becomes smaller than a few wavelengths of the light.
In this case, the distance between two scattering events is so short that (1)
each scattered photon does not reach the far field limit before being re-
scattered and (2) the (differential) scattering cross section of each and
every scatterer is affected by multiply scattered photons returning to it. A
first attempt to describe these near-field effects was recently proposed by
Rezvani Naraghi  _et al._ (2015), but takes into account only the first point.
In this paper, we propose a different light transport model that takes into
account both effects by considering that each scatterer in a densely packed
random assembly is embedded in a properly estimated effective medium 111The
importance of an appropriate description of the effective refractive index to
describe light transport in randomly packed colloids was already noted by
Reufer  _et al._ (2007).. We show that the use of the energy-density coherent
potential approximation (ECPA, Busch and Soukoulis (1995); *Busch1996) for the
calculation of the effective refractive index provides a proper description of
light transport in densely packed multiple scattering samples. The predictions
of this model are tested experimentally on samples having resonant transport
properties, namely, monodisperse photonic glasses, which were first introduced
and studied by García and coworkers García  _et al._ (2007); Sapienza  _et
al._ (2007); García  _et al._ (2008). These materials are made of randomly
assembled monodisperse Mie scatterers (which are isolated dielectric spheres
of the same radius $r$). A common idea is that single scattering Mie
resonances (see Fig. 1(a) or Ref. Bohren and Huffman, 1998) have a signature
in the transport mean free path $\ell^{\star}$ of randomly assembled Mie
scatterers. These resonances were already observed, but up to now, all
attempts to describe them more than qualitatively failed Sapienza  _et al._
(2007); Reufer  _et al._ (2007); Montesdeoca  _et al._ (2016). We argue here
that a proper description of resonant near-field effects is the key to a
quantitative description of them.

In this paper, we present a fit parameter free model that is able to predict
for the first time the correct positions ($r/\lambda_{0}$) of the resonance
peaks of the scattering strength ($\lambda_{0}/\ell^{\star}$) and describes
their order of magnitude. We tested this model against ab initio numerical
simulations, earlier experimental data obtained from transmission experiments
García  _et al._ (2008), and new experimental results from backscattering
experiments on specially synthesized polystyrene colloidal glasses Chen  _et
al._ (2017).

Our description of the light transport in photonic glasses starts with the
derivation of the transport mean free path Wolf  _et al._ (1988),

$\ell^{\star}=\frac{\ell_{\mathrm{s}}}{1-\langle\cos\theta\rangle}\overset{\text{sphere}}{=}\frac{1}{1-\langle\cos\theta\rangle}\frac{4\pi
r^{3}}{3f\sigma_{\mathrm{s}}}.$ (1)

The second equality holds in the case of an assembly of monodisperse Mie
scatterers. $\ell_{\mathrm{s}}$ is the scattering mean free path (the distance
between two scattering events), $f$ is the filling fraction of the scatterers,
and $\theta$ is the scattering angle. Both the scattering cross section
$\sigma_{\mathrm{s}}$ and the anisotropy factor $\langle\cos\theta\rangle$ of
each scattering event in the multiple scattering regime can be expressed in
terms of the form factor $F(\theta)$ (scattering properties of the single
sphere) and the structure factor $S(\theta)$ (collective scattering of the
sample) Fraden and Maret (1990); van de Hulst (1957)

$\sigma_{\mathrm{s}}=\frac{\pi}{k^{2}}\int_{0}^{\pi}F(\theta)S(\theta)\sin\theta\,\mathrm{d}\theta,$
(2) $\langle\cos\theta\rangle=\dfrac{\int_{0}^{\pi}\cos\theta
F(\theta)S(\theta)\sin\theta\,\mathrm{d}\theta}{\int_{0}^{\pi}F(\theta)S(\theta)\sin\theta\,\mathrm{d}\theta}.$
(3)

The intensity form factor $F(\theta)$ of a sphere is calculated by the Mie-
theory Bohren and Huffman (1998). Since the samples we focus on are assemblies
of randomly packed spheres, we use the hard sphere Percus-Yevick structure
factor Percus and Yevick (1958) $S(q)$ with $q=2k\sin\theta/2$ for $S(\theta)$
($k=2\pi/\lambda$ with $\lambda$ the wavelength of the light in the
surrounding medium). Using this description, Fraden and Maret (1990) were able
to describe the effect of short-range correlations of the scatterers positions
on multiple scattering of light in polystyrene spherical colloids ($r=230$ nm)
suspended in water with filling fractions up to 45% (refractive indices
$n_{\mathrm{PS}}=1.60$ in $n_{\mathrm{H_{2}O}}=1.33$).

However, when the average distance between nearest colloids is of the order of
the light wavelength, near-field effects start to play a role in the transport
properties Rezvani Naraghi  _et al._ (2015). These near-field effects are not
caught in $F(\theta)$ when calculated with the bulk refractive index of the
surrounding medium. Moreover, $S(q)$ depends on the surrounding refractive
index via $\lambda$ as well. As we show below, the presence of other
scatterers in the direct vicinity can be taken into account by defining an
effective refractive index $n_{\mathrm{eff}}$ for the surrounding medium.
$n_{\mathrm{eff}}$, which can first be roughly estimated as the volume average
of the local refractive index, lowers the refractive index difference between
each scatterer and the medium, and therefore lowers its scattering strength.
This effect was not considered in Ref. Fraden and Maret, 1990, but does not
alter much that analysis because of the relatively low refractive index
contrast between polystyrene and water. We show here that a proper estimation
of $n_{\mathrm{eff}}$ is necessary when dealing with higher refractive index
contrasts (like polystyrene in air, $n_{0}=1$) 222Let us just mention here
that we are not aware of any direct measurement of this quantity..

To estimate the effective refractive index of our photonic glasses, we use the
ECPA by Busch and Soukoulis (1995) 333A comparison between the ECPA and other
ways of defining $n_{\mathrm{eff}}$ van Albada  _et al._ (1991); Soukoulis
_et al._ (1994, 1989); Economou and Soukoulis (1989) is given in Appendix A..
Since the positions of the scatterers are random (with the constraint that the
spheres can not overlap), each one will be coated, on average, by a shell
having the refractive index of the embedding matrix. This _core-shell_
particle itself is embedded in a medium having the effective refractive index
(see inset of Fig. 1(b)). By definition, the energy density inside an
effective medium should be uniform: we find $n_{\mathrm{ECPA}}$ using an
iterative process such that the energy density inside the core-shell particle
matches the energy density inside the same volume, but having the uniform
refractive index $n_{\mathrm{ECPA}}$. The result of this iterative process is
plotted in Fig. 1(b) as a function of the size parameter and of the scatterer
density for polystyrene colloids embedded in air.

Figure 1: (a) Scattering efficiency $Q_{\mathrm{sca}}$ of an isolated Mie
scatterer having a refractive index of $n_{\mathrm{PS}}=1.6$ in air. The
shaded regions highlight the region with negative curvature and the dashed
white lines are a guide to the eye. – (b) Color map of $n_{\mathrm{ECPA}}$
calculated for a refractive index of the particle $n_{\mathrm{PS}}$ surrounded
by air. Inset: principle of the ECPA (see text) – (c) Color map of the
scattering strength $\lambda_{0}/\ell^{\star}$ calculated with Eq. (1) using
$n_{\mathrm{ECPA}}$. In both cases, the quantities are shown as a function of
both the size parameter $r/\lambda_{0}$ and the filling fraction $f$.

Clear resonances are seen as a function of the size parameter. We can
furthermore see an increase of the refractive index with increasing filling
fraction. We note, finally, that the largest $n_{\mathrm{ECPA}}$ corresponds
approximately to the index of the scattering particles
($n_{\mathrm{PS}}=1.6$).

Replacing the wave vector $k$ by $k_{\mathrm{eff}}=2\pi
n_{\mathrm{eff}}/\lambda_{0}$ in the calculation of $\sigma_{\mathrm{s}}$ and
$\langle\cos\theta\rangle$, we can calculate the scattering strength
$\lambda_{0}/\ell^{\star}$ using Eq. (1) as a function of both $r/\lambda_{0}$
and $f$ (see Fig. 1(c)). Again resonances can be seen as a function of the
size parameter. Note that the peaks of $\lambda_{0}/\ell^{\star}$ are placed
at dips of $n_{\mathrm{ECPA}}$. This can be understood in terms of partial
index matching between the scatterers and the effective surrounding medium.
Fig. 2(a) shows the scattering strength for a given filling fraction ($f=0.5$,
red solid curve) and compares it to the one expected taking a given
polydispersity of the spheres into account (dash-dotted red curve).

Figure 2: (a) Scattering strength $\lambda_{0}/\ell^{\star}$ calculated with
different models for a filling fraction $f=0.5$ (red curves: using
$n_{\mathrm{ECPA}}$ (solid: monodisperse, dot dashed: 5% polydispersity; the
shaded regions highlight negative curvature of the solid red curve); dashed
curves: simple model for $n_{\mathrm{eff}}$ (see text); dotted curves:
calculated with the differential cross section of the $N$ cluster). The single
sphere scattering efficiency (also plotted in Fig. 1(a) is shown for
comparison in gray (arbitrary units). – (b) Averaged
$\lambda_{0}/\ell_{\mathrm{s}}$ calculated with the MSTM code over five
realizations of numerical photonic glasses, each one having $f=0.5$.

As expected, the resonances are somewhat smeared out.

#### Comparison with other models:

To emphasize our claim that a proper choice of $n_{\mathrm{eff}}$ is necessary
for the calculation of the transport properties of dense high index photonic
glasses, we now show the difference between the scattering strength calculated
using $n_{\mathrm{ECPA}}$ (red solid curve in Fig. 2) and the one calculated
using more simple models ($n_{\mathrm{vol}}=f\cdot n_{\mathrm{PS}}+(1-f)\cdot
n_{0}$ (violet) or $n_{\mathrm{MG}}$ for the Maxwell Garnett effective
refractive index Garnett (1904) (green)). The simple models lead to no
resonant behavior in $n$ but the order of magnitude of the scattering strength
agrees more with the $n_{\mathrm{ECPA}}$ predictions than the curve calculated
without using any effective refractive index (blue dashed curve).

An other way to take the presence of other scatterers in the direct vicinity
of each very scatterer into account for its scattering properties is to
calculate directly $F(\theta)$ for a scatterer surrounded by other scatterers
(as suggested by García  _et al._ (2008)). We therefore calculated numerically
using the Multiple Sphere T Matrix (MSTM) code Mackowski and Mishchenko (2011)
the average over all possible orientations of the differential scattering
cross section $F_{N}(\theta,d_{\mathrm{NN}}(f))$ of clusters of $N=2$, 7, or
13 particles made of one sphere having $N-1$ neighboring spheres placed at the
average nearest-neighbor distance in the glass $d_{\mathrm{NN}}$ (which is a
function of $f$ Torquato (1995)). We then replaced $F(\theta)$ by
$F_{N}(\theta,d_{\mathrm{NN}}(f))/N$ in Eqs. (2) and (3) _without_ using an
effective refractive index. By doing this, we take into account near-field
effects in the scattering properties in a similar way as in Ref. Rezvani
Naraghi  _et al._ , 2015. As shown by the dotted curves in Fig. 2(a), the
scattering strength decreases with the number of nearest neighbors indicating
stronger near-field effects, but is still different than the one predicted
using $n_{\mathrm{ECPA}}$.

#### Numerical test:

To test the model using $n_{\mathrm{ECPA}}$, we simulate the transmission
trough thin slabs (6.25 times the radius of the spheres) of monodisperse
spheres using the MSTM code Mackowski and Mishchenko (2011). The glasses have
a filling fraction of 50%, and were numerically synthesized by using a force-
biased algorithm Baranau and Tallarek (2014). In our simulations, we
calculated the transmission of the non-scattered field and extracted the
_scattering_ mean free path from this quantity 444See Appendix B.. As
expected, the so-extracted $\ell_{\mathrm{s}}$ are smaller than the
$\ell^{\star}$ predicted by our model. Moreover, one can recognize the same
resonant behavior between $\ell_{\mathrm{s}}$ and $\ell^{\star}$, at least for
$r/\lambda_{0}\lesssim 0.7$. For larger ones, the resonances seems to be
smeared out (in this range, the overall transmissions fall in the $10^{-3}$
range, making it difficult to extract the coherent beam without any scattered
contribution).

#### Experimental test:

To test the model experimentally, we measure the transport mean free path
$\ell^{\star}$ of photonic glasses for different size parameters
$r/\lambda_{0}$ by analyzing the shape of their coherent backscattering cone
(CBC) Wolf and Maret (1985); Albada and Lagendijk (1985). The width of the CBC
is inversely proportional to $k\ell^{\star}$. In this method a parallel light
beam illuminates the multiple scattering sample via a beamsplitter (see Fig.
3), and the backscattered light is imaged in the Fourier space on a CCD camera
(Apogee Alta U4000) placed in the focal plane of a convex lens
($f^{\prime}=200\,\mathrm{mm}$).

Figure 3: Scheme of the CBC setup used to measure $\ell^{\star}$ of photonic
glasses: Illumination can be tuned from 450 nm to 700 nm via a white LED
coupled in a monochromator. The output illuminates the sample via a 50:50
beamsplitter and a circular polarizer is used to filter single scattered
light. The reflected light is imaged in the focal plane of a $200$ mm lens on
a CCD camera. A typical image of the cone of a polystyrene photonic glass is
shown on the left of the CCD. An image of a macroscopic polystyrene sample is
shown in the lower left with a SEM image of the same sample to its left.

A circular polarizer in front of the sample filters single backscattered light
555In single scattering events the helicity is flipped and filtered while
multiple scattering randomizes the polarization.. This setup is able to
measure angles up to $3.0^{\circ}$. To probe $\ell^{\star}$ for different
$r/\lambda_{0}$ ratios, we used a white LED 666Lumileds LUXEON CoB 109 5000 K
coupled in a monochromator 777Acton SP-2150i, Princeton Instruments as light
source to tune $\lambda_{0}$ between 450 and 700 nm with a wavelength width of
$\Delta\lambda_{0}\approx 5\,\mathrm{nm}$. Since this is an incoherent source,
there is no need to average over speckles to measure the CBC 888This light
source has an slightly divergent beam which needs to be accounted for in the
analysis, see Appendix C.. Our samples are free standing photonic glasses
having a thickness $L\approx 1\,\mathrm{mm}$ and a diameter $d\approx
1\,\mathrm{cm}$. The samples were prepared from monodisperse polystyrene
colloids in solution via ultracentrifugation and with addition of a controlled
amount of salt Chen  _et al._ (2017) to avoid crystallization and prevent
optical shortcuts. A small amount of polyacrylamide holds the particles
together to allow freestanding photonic glasses. Therefore, from the
measurement of the mass and of the volume of the samples, we can only estimate
an upper limit for the filling fractions of the particles in the glasses.
Their filling fraction might be slightly lower than 0.55. With our light
source and by varying the particles size from $r=125$ to $r=335\,\mathrm{nm}$
999The radii and their errors are obtained from the analysis of SEM images (as
the one shown in Fig. 3 for the sample with the largest particles). All
samples have polydispersities lower than $5\%$. we are able to cover a wide
range of $r/\lambda_{0}$ from 0.17 to 0.70.

The image in Fig. 3 shows the CBC of a polystyrene photonic glass. For each
measurement we averaged about five images (exposure time
$t_{\mathrm{exp}}=3\,\mathrm{s}$) before performing a radial average. The
obtained data were finally fitted with the standard CBC formula to obtain
$\ell^{\star}$ Akkermans and Montambaux (2007). The diffusion constant $D$ and
the absorption time $\tau_{\mathrm{a}}$ were measured separately in time of
flight experiments to extract the absorption length
$L_{\mathrm{a}}=\sqrt{D\tau_{\mathrm{a}}}$ which is used in the CBC fit.
$L_{\mathrm{a}}$ was typically on the order of the sample thickness
$L\gg\ell^{\star}$, meaning that absorption was very low. A wavelength scan
for five different photonic glasses is shown in Fig. 4 and reveals the first
direct observation of strong resonances of $\lambda_{0}/\ell^{\star}$ in the
visible.

Figure 4:  Measured scattering strength of five monodisperse polystyrene
photonic glasses from $r=125$ to $335\,\mathrm{nm}$ compared to a chalk sample
(irregular particle shape, for this sample we arbitrarily set $r=335$ nm). The
errors in $\ell^{\star}$ are estimated from the fits to be $\pm
0.1\,\mu\mathrm{m}$. The data are compared with the predictions of the
scattering model (Eq. (1) using $n_{\mathrm{ECPA}}$) for different filling
fractions (solid lines) and with the transmission data extracted in García
_et al._ (2008) ($r=610$ nm).

For comparison we measured a randomly shaped and highly polydisperse sample of
chalk powder (black diamonds) and, as expected, no resonant behavior is
observed in this case. The experimental data show that the positions of the
resonances are very well predicted when the transport properties are
calculated using the ECPA effective refractive index (see the comparison
between the scattering strengths predicted by different models in Fig. 2). The
amplitudes of the resonances are slightly smeared out compared to the model.
This can partially be understood because of the residual polydispersity of the
particles and of the wavelength width of the source (see Fig. 2). To further
test our model we plotted the data of Fig. 19 of García  _et al._ (2008) (red
diamonds). They measured $\ell^{\star}$ via diffuse transmission measurements
in the infrared with polystyrene particles of $r=610\,\mathrm{nm}$. Our model
recovers the position of their resonances very well. The lower amplitude of
the measured resonances compared with the model could be explained by the
polydispersity of their sample, but no such information is given in the paper.
All these six different data sets collapse on the same curve and cover the
first nine Mie resonances versus $r/\lambda_{0}$, and this without any fit
parameter.

In this paper we highlight the importance of an appropriate model for
$n_{\mathrm{eff}}$ in the description of the transport properties of dense
photonic materials where near-field coupling influences scattering behavior.
The shell introduced in the definition of $n_{\mathrm{ECPA}}$ is the key to
take near-field effects into account. It has a thickness which is related to
the average particle distance in the glass, and couples electromagnetically
each scatterer with the surrounding medium. Unlike all other models used so
far, the ECPA model predicts very well the resonant behavior in photonic
glasses without any fit parameter.

The resonant behavior in polymeric photonic glasses leads to strong scattering
at certain $r/\lambda_{0}$ values. Increasing the refractive index contrast by
using higher index materials such as TiO2 spheres in air will lead to much
stronger resonances. While the preparation of freestanding photonic glasses
made of such high index colloids with the right $r/\lambda_{0}$ values may
still be a technical issue, it would possibly give access to a new scattering
regime where signatures of three-dimensional light localization are expected.

###### Acknowledgements.

The authors would like to thank Ralf Tweer for providing his ECPA-code and
Herbert Kaiser for useful input on the numerical part of this paper. We
further thank the Center for Applied Photonics (CAP) at the Universität
Konstanz and the Schweizerischer Nationalfonds (SNF) for financial support.
The numerical simulations were done on the Scientific Compute Cluster at the
Universität Konstanz. GJA acknowledges support from the Zukunftskolleg
(Universität Konstanz) for an Independent Research Start-up Grant.

##  Appendix A Effective refractive index

Beside the need of an effective index theory that takes into account the
resonant behavior of Mie scatterers (like the subject of this paper),
different effective medium theories were developed in the quest for an
appropriate description of transport properties in the high concentration
regime in multiple scattering samples.

The energy transport velocity $v_{\mathrm{E}}$ of a wave propagating in a
random medium is lowered with respect to the speed of light in vacuum,
$v_{\mathrm{E}}=\frac{c}{n_{\mathrm{eff}}}$, $n_{\mathrm{eff}}$ being the
effective refractive index of the complex material. For very small
concentrations of scattering particles $n_{\mathrm{eff}}$ tends towards the
refractive index of the matrix medium $n_{\mathrm{m}}$. If the sample is
completely filled by the particle medium $n_{\mathrm{eff}}$ converges to
$n_{\mathrm{p}}$. The easiest way to define an effective refractive index
which captures this behavior is

$n_{\mathrm{eff}}=fn_{\mathrm{p}}+(1-f)n_{\mathrm{m}}.$ (4)

This is only valid for extreme situations of very low or very high filling
fractions $f$.

Another way to calculate the effective refractive index is the Maxwell Garnett
mixing formula Garnett (1904)

$\displaystyle n_{\mathrm{MG}}=n_{m}\sqrt{1+\frac{3fK}{1-fK}}$ (5)

with
$K=\frac{n_{\mathrm{p}}^{2}-n_{\mathrm{m}}^{2}}{n_{\mathrm{p}}^{2}+2n_{\mathrm{m}}^{2}}$.
The latter model is more physical than Eq. (4) because it assumes that the
polarizabilities are additive, not the refractive indexes. Nevertheless, even
if the derivation of Eq. (5) assumes the scatterers to be spherical, it
completely neglects resonant scattering (the fact that energy can be stored in
Mie scatterers for certain radius over wavelength ratios).

By using a Bethe-Salpeter equation, van Albada  _et al._ (1991) were able to
describe the effect of these Mie resonances in the regime of low $f$, and
their theory was consistent with experiments at $f=0.36$. This showed that
experimentally obtained low values of the diffusion constant
$D=\frac{v_{\mathrm{E}}\ell^{\star}}{3}$ were related to low values of
$v_{\mathrm{E}}$ due to resonant Mie scattering and a corresponding energy
storage process, and not to low values of $\ell^{\star}$ which would signify
localization Busch and Soukoulis (1995); *Busch1996. For higher filling
fractions (like in the case studied in this paper) experimental data are much
better predicted by advanced versions of the so called coherent potential
approximation (CPA) Soukoulis  _et al._ (1994).

###  A.1 Coherent Potential Approximation

In the “genuine” CPA Soukoulis  _et al._ (1989); *Economou1989 a medium of two
lossless materials consisting of spheres with radius $r$, refractive index
$n_{\mathrm{p}}$, and a volume fraction $f$ in a host material with refractive
index $n_{\mathrm{m}}$ is considered. Each point in the medium will be either
in a region with $n_{\mathrm{p}}$ with a probability $f$, or in a region with
$n_{\mathrm{m}}$ with a probability $1-f$. The medium is modeled by spheres of
radius $r$ having a refractive index given by the aforementioned
probabilities. The effective refractive index is then found such that the
averaged forward scattering amplitude vanishes,

$\displaystyle fF_{\mathrm{p}}(0)+(1-f)F_{\mathrm{m}}(0)=0,$ (6)

with $F_{\mathrm{p}}$ (resp. $F_{\mathrm{m}}$) being the differential
scattering cross-section of a sphere of refractive index $n_{\mathrm{p}}$
(resp. $n_{\mathrm{m}}$) embedded in a medium having the effective refractive
index.

This approach neglects topological and geometrical differences between the
scattering spheres and the host material, e.g., for high $f$ the host spheres
are not only less probable but also have to have smaller radii. The real
random system would be much better estimated by a mixture of coated sphere
(where the scatterer with $n_{\mathrm{p}}$ is coated by a spherical region of
host material $n_{\mathrm{m}}$) and of spheres (of refractive index
$n_{\mathrm{m}}$). Equation (6) (this time with $F_{\mathrm{p}}$ the
differential scattering cross-section of the coated sphere) is then solved in
the same way as in the classical CPA to obtain the effective medium refractive
index $n_{\mathrm{eff}}$ self-consistently. In this approach the coating
thickness varies with $f$. Due to the condition that the spheres should not
overlap, the distribution of spacings between neighboring spheres has a peak
at $r_{\mathrm{c}}>r$ with $r_{\mathrm{c}}=\frac{r}{f^{1/3}}$. This advanced
version of the CPA has the advantage of taking into account short range order
and thus fits to experimental data quite well in the high $f$ regime Soukoulis
_et al._ (1994). Nevertheless, the effective refractive index calculated by
the coated CPA is very close to the one calculated with the Maxwell-Garnett
theory (see Fig. 4 of Ref. Soukoulis  _et al._ , 1994), whereas Fig. 2(a)
shows that the Maxwell-Garnett theory does not fit our data.

###  A.2 Energy-Density Coherent Potential Approximation

In the classical CPA the energy density is homogeneous by construction. This
can be violated in the coated CPA approach because of the coated sphere as the
basic scattering unit. For low $f$ (large coatings) this leads to unphysical
behavior near the single sphere Mie resonances, e.g., refractive indices
smaller than 1 such that $v_{\mathrm{E}}>c$. Therefore a new CPA approach was
developed by Busch and Soukoulis (1995), the so-called energy-density coherent
potential approximation. Here the termination criterion for the iterative
determination of the effective refractive index $n_{\mathrm{eff}}$ is that a
homogeneous energy density $\rho_{\mathrm{E}}$ on scales larger than the basic
scattering unit is reached. This is schematically shown in the inset of Fig.
1(b).

The criterion of a constant energy density in the case of a plane wave hitting
a coated sphere embedded in the effective medium versus the case where the
same volume is filled by the effective medium only is quantitatively expressed
in the self-consistent equation

$\int_{0}^{r_{\mathrm{c}}}\mathrm{d}^{3}R\rho^{(1)}_{\mathrm{E}}(\vec{R})=\int_{0}^{r_{\mathrm{c}}}\mathrm{d}^{3}R\rho^{(2)}_{\mathrm{E}}(\vec{R})\
,$ (7)

where $\rho^{(1)}_{\mathrm{E}}(\vec{R})$ and
$\rho^{(2)}_{\mathrm{E}}(\vec{R})$ are the energy densities in the coated
sphere and in the same volume filled with a medium having a refractive index
$n_{\mathrm{eff}}$ respectively. The energy density of a electromagnetic wave
can be expressed as:

$\rho_{\mathrm{E}}(\vec{R})=\frac{1}{2}\left[\varepsilon(\vec{R})|\vec{E}(\vec{R})|^{2}+\mu|\vec{H}(\vec{R})|^{2}\right]\
.$ (8)

where $\vec{E}$ and $\vec{H}$ are the electric and magnetic fields,
$\varepsilon$ is the dielectric constant, and $\mu$ is the magnetic
permeability (the latter is assumed to be the same in both materials). With
these equations the effective refractive index can be determined for all
frequencies guaranteeing a homogeneous energy density on scales larger than
the scattering unit.

Busch and Soukoulis (1995) claim that their model, which takes into account
multiple scattering effects in a mean-field sense, can be used for scalar
(acoustic), vector (electromagnetic), and tensor (elastic) waves and is valid
for all densities of scatterers. They also test it on earlier experimental
data. More interesting for the present paper, this model captures the effects
of resonant near-field effects in the multiple scattering regime.

##  Appendix B Numerical simulations

The transmission calculations of Fig. 2(b) were obtained by using the Multiple
Sphere T Matrix code Mackowski and Mishchenko (2011). The geometry is shown in
Fig. 5: we created cylindrical slabs having a diameter of 10000 nm, a
thickness of $L=1000$ nm containing about 2270 particles with a radius $r=160$
nm (filling fraction of 50%).

Figure 5: MSTM geometry: a Gaussian beam is focused on one side of the sample,
and the intensity is calculated in the far field on the red spot (see text for
details). The figure is at scale.

The five samples were created by using the force-biased algorithm code of
Vasili Baranau Baranau and Tallarek (2014).

A Gaussian beam (wavelength $\lambda_{0}$ between 160 and 1600 nm) was focused
on one side of the slab (waist $w_{0}=1500$ nm), and $I=\left<E^{2}\right>$
—the integral of the intensity on a circle of radius101010The influence of the
target size is shown in Fig. 6. 100 nm placed at a distance of 16000 nm on the
other side of the sample (ten times the largest wavelength to get rid of any
evanescent wave present close to the spheres)— was calculated111111Five
different glass configurations were calculated on the Scientific Compute
Cluster (Universität Konstanz) using between 10 and 20 processors for
computing time ranging from a few hours to a few days depending on the
wavelength.. We did the same calculation without any scatterers
($I_{0}=\left<{E_{0}}^{2}\right>$) and the corresponding transmission values
are shown in Fig. 6 for different target sizes. These samples are optically
thin ($\lambda_{0}/L$ ranges from 0.16 to 1.6). We therefore expect most of
the photons to be scattered only once. In this regime, the coherent part of
the beam $I_{\mathrm{c}}$ (i.e. the part of the wave which is not scattered)
is attenuated exponentially

$\displaystyle
I_{\mathrm{c}}=I_{0}\cdot\exp{\left(-\frac{L}{\ell_{\mathrm{s}}}\right)},$ (9)

where $I_{0}$ is the incident intensity and $\ell_{\mathrm{s}}$ is the
_scattering_ mean free path.

Figure 6: Average over five different slabs of the MSTM calculated
transmission as a function of the size parameter: influence of the target size
on the transmission values. The target is either circular (violet, green, blue
and yellow lower curves) or square (red upper curve, square area $10640^{2}$
nm2). Away from the focus, the width of a Gaussian beam depends on the
wavelength: the yellow dashed curve was calculated by taking the radius of the
target such that $I_{0}(r)>I_{0}(0)/2$. All error bars correspond to the
standard deviation.

By integrating the intensity on a small surface at such a large distance on
such a thin sample, the major part of it corresponds to the coherent
intensity, therefore $I\simeq I_{\mathrm{c}}$. This allows us to calculate
$\ell_{\mathrm{s}}$ using Eq. (9).

##  Appendix C Data analysis with divergence

|
---|---
Figure 7: Sketch of the effect of the divergence on a focused beam and
measured image of the divergent spot.

Our light source (LED + monochromator) has a divergent beam compared to a
coherent laser beam. This needs to be accounted for in the analysis of the
coherent backscattering data. The backscattered light from the sample diverges
slightly before being captured by the lens and imaged in the Fourier space.
This leads to a smeared out (in angles) cone measurement. The divergence of
the parallelized beam from the LED compared to a perfectly parallel beam can
be measured by replacing the sample by a mirror. For a parallel beam this
leads to a perfectly focused spot on the CCD if the setup is well adjusted.
The spot imaged from a divergent beam is wider and reflects the angular
divergence of the light source (see Fig. 7). The spot has a radius of
$r\approx 1.4\leavevmode\nobreak\ \mathrm{mm}$. For the analysis a binary
image was generated with a circular spot of the same size centered around the
center of the measured CBC. This calibration image was then convoluted with
the CBC formula Akkermans and Montambaux (2007) as follows: a two-dimensional
image was generated from the CBC formula and transformed with a discrete
cosine transform (DCT, a similar function as a discrete Fourier transform but
having different boundary conditions). The calibration image was also
transformed by a DCT and both functions are multiplied in Fourier space. The
result is then transformed back by an inverse DCT. This function is finally
used to fit the measured image of the CBC.

## References

  * Abrahams  _et al._ (1979) E. Abrahams, P. W. Anderson, D. C. Licciardello, and T. V. Ramakrishnan, Phys. Rev. Lett. 42, 673 (1979).
  * Bayer and Niederdränk (1993) G. Bayer and T. Niederdränk, Phys. Rev. Lett. 70, 3884 (1993).
  * Larose  _et al._ (2004) E. Larose, L. Margerin, B. A. van Tiggelen, and M. Campillo, Phys. Rev. Lett. 93, 048501 (2004).
  * Jendrzejewski  _et al._ (2012) F. Jendrzejewski, K. Müller, J. Richard, A. Date, T. Plisson, P. Bouyer, A. Aspect, and V. Josse, Phys. Rev. Lett. 109, 195302 (2012).
  * Wolf and Maret (1985) P.-E. Wolf and G. Maret, Phys. Rev. Lett. 55, 2696 (1985).
  * Albada and Lagendijk (1985) M. P. V. Albada and A. Lagendijk, Phys. Rev. Lett. 55, 2692 (1985).
  * Akkermans and Montambaux (2007) E. Akkermans and G. Montambaux,  _Mesoscopic Physics of Electrons and Photons_ (Cambridge University Press, 2007).
  * Lagendijk  _et al._ (2009) A. Lagendijk, B. van Tiggelen, and D. S. Wiersma, Physics Today 62, 24 (2009).
  * Fraden and Maret (1990) S. Fraden and G. Maret, Phys. Rev. Lett. 65, 512 (1990).
  * Rezvani Naraghi  _et al._ (2015) R. Rezvani Naraghi, S. Sukhov, J. J. Sáenz, and A. Dogariu, Phys. Rev. Lett. 115, 203903 (2015).
  * Reufer  _et al._ (2007) M. Reufer, L. F. Rojas-Ochoa, S. Eiden, J. J. Sáenz, and F. Scheffold, Applied Physics Letters 91, 171904 (2007).
  * Busch and Soukoulis (1995) K. Busch and C. M. Soukoulis, Phys. Rev. Lett. 75, 3442 (1995).
  * Busch and Soukoulis (1996) K. Busch and C. M. Soukoulis, Phys. Rev. B 54, 893 (1996).
  * García  _et al._ (2007) P. García, R. Sapienza, Á. Blanco, and C. López, Advanced Materials 19, 2597 (2007).
  * Sapienza  _et al._ (2007) R. Sapienza, P. D. García, J. Bertolotti, M. D. Martín, Á. Blanco, L. Viña, C. López, and D. S. Wiersma, Phys. Rev. Lett. 99, 233902 (2007).
  * García  _et al._ (2008) P. D. García, R. Sapienza, J. Bertolotti, M. D. Martín, A. Blanco, A. Altube, L. Viña, D. S. Wiersma, and C. López, Phys. Rev. A 78, 023823 (2008).
  * Bohren and Huffman (1998) C. F. Bohren and D. R. Huffman,  _Absorption and Scattering of Light by Small Particles_ (Wiley, New York, 1998).
  * Montesdeoca  _et al._ (2016) D. Montesdeoca, F. Bayat, A. Espinha, Á. Blanco, C. Pecharromán, and C. López, Particle & Particle Systems Characterization 33, 871 (2016).
  * Chen  _et al._ (2017) M. Chen, D. Fischli, L. Schertel, G. J. Aubry, B. Häusele, S. Polarz, G. Maret, and H. Cölfen, Small 13, 1701392 (2017).
  * Wolf  _et al._ (1988) P. Wolf, G. Maret, E. Akkermans, and R. Maynard, Journal de Physique 49, 63 (1988).
  * van de Hulst (1957) H. van de Hulst,  _Light Scattering by Small Particles_, Dover Books on Physics Series (Dover Publications, New York, 1957).
  * Percus and Yevick (1958) J. K. Percus and G. J. Yevick, Phys. Rev. 110, 1 (1958).
  * van Albada  _et al._ (1991) M. P. van Albada, B. A. van Tiggelen, A. Lagendijk, and A. Tip, Phys. Rev. Lett. 66, 3132 (1991).
  * Soukoulis  _et al._ (1994) C. M. Soukoulis, S. Datta, and E. N. Economou, Phys. Rev. B 49, 3800 (1994).
  * Soukoulis  _et al._ (1989) C. M. Soukoulis, E. N. Economou, G. S. Grest, and M. H. Cohen, Phys. Rev. Lett. 62, 575 (1989).
  * Economou and Soukoulis (1989) E. N. Economou and C. M. Soukoulis, Phys. Rev. B 40, 7977 (1989).
  * Garnett (1904) J. C. M. Garnett, Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 203, 385 (1904).
  * Mackowski and Mishchenko (2011) D. Mackowski and M. Mishchenko, Journal of Quantitative Spectroscopy and Radiative Transfer 112, 2182 (2011).
  * Torquato (1995) S. Torquato, Phys. Rev. Lett. 74, 2156 (1995).
  * Baranau and Tallarek (2014) V. Baranau and U. Tallarek, Soft Matter 10, 3826 (2014), code available at https://doi.org/10.5281/zenodo.580324.