The Tail of Late-forming Dwarf Galaxies in $\Lambda$CDM

We use a robust analytical model together with a high-resolution hydrodynamical cosmological simulation to demonstrate that in a $\Lambda$ cold dark matter ($\Lambda$CDM) universe, a small fraction of dwarf galaxies inhabiting dark matter (DM) halos in the mass range $3\times 10^{9} \lesssim M_{200} / M_{\odot} \lesssim 10^{10}$ form unusually late ($z<3$) compared to the bulk population of galaxies. These galaxies originate from the interplay between the stochastic growth of DM halos and the existence of a time-dependent DM halo mass below which galaxies do not form. The formation epoch of the simulated late-forming galaxies traces remarkably well the time when their host DM halos first exceeded a nontrivial (but well-understood) time-dependent critical mass, thus making late-forming dwarfs attractive cosmological probes with constraining power over the past growth history of their host halos. The agreement between our model and the simulation results demonstrates that the population of simulated late-forming dwarfs is a robust cosmological outcome and largely independent of the specific galaxy formation model included in the simulations provided: (1) the universe underwent cosmic reionization before $z_{\rm re} \sim 8$; (2) star formation proceeds in gas that self-gravitates; and (3) galaxy formation is largely restricted to atomic-cooling halos before $z_{\rm re}$. The scarcity of massive late-forming dwarfs expected in $\Lambda$CDM implies that the great majority of bright, metal-poor, and actively star-forming dwarfs observed in our local universe--the most obvious candidates for these late-forming galaxies--cannot be undergoing their formation for the first time at the present day in a $\Lambda$CDM universe.


INTRODUCTION
Within the Lambda cold dark matter (ΛCDM) cosmological model, galaxies form from gas that collapses in the center of gravitationally bound dark matter (DM) halos (White & Rees 1978). Galaxies do not form, however, in every halo. In the absence of external heating sources, galaxy formation is restricted to the so-called atomic-cooling (AC) halos, i.e., halos that shock heat the gas to a temperature T 10 4 K, above which radiative cooling becomes efficient. Under the presence of external energetic sources that suppress atomic cooling and heat the gas (e.g. Efstathiou 1992), such as the ultraviolet background (UVB) radiation field that keeps the universe (re)ionized, galaxy formation proceeds in halos that exceed a time-dependent critical mass above which the gas becomes self-gravitating. The value of this critical mass has been estimated in the past using either idealized Jeans arguments or numerical simulations (e.g. Rees 1986; Thoul & Weinberg 1996;Quinn et al. 1996), 1 leading to the under-standing that for galaxies to form after cosmic reionization (CR), their host halos must exceed the AC limit by a factor of a few. The AC limit is, therefore, useful to determine which halos host galaxies at high redshift, before the universe underwent CR (e.g. Oh & Haiman 2002;Xu et al. 2013;Wise et al. 2014;Xu et al. 2016), whereas the critical mass describes the onset of galaxy formation after CR. The onset of galaxy formation is thus deeply linked to the growth of DM halos, as only those halos that exceed the critical mass can host a luminous galaxy in their center. These ideas have been shown to agree qualitatively with results of hydrodynamical cosmological simulations (e.g. Okamoto & Frenk 2009;Sawala et al. 2016;Benítez-Llambay et al. 2015Fitts et al. 2017;Macciò et al. 2017, and references therein), and are of fundamental importance to explain the scarcity of observed luminous satellites compared to results of collisionless simulations (Klypin et al. 1999;Bullock et al. 2000). However, it is clear that neither halos have uniform density nor the interstellar medium is isothermal, so the successful characterization of the critical mass necessarily requires more advanced modeling. Benitez-Llambay & Frenk (2020) (hereafter BLF) have recently derived the critical mass for the onset of galaxy formation considering the nonlinear collapse of DM halos and the distinctive temperature-density relation of the intergalactic medium, removing the freedom inherent to Jeans arguments. Their critical mass (hereafter BLF mass), which differs from that arising from idealized Jeans modeling, is remarkably accurate compared to results of nonlinear hydrodynamical cosmological simulations.
The existence of a critical mass for galaxy formation, coupled with the growth history of DM halos, has some interesting consequences. First, galaxies must form stochastically in halos of present-day mass under the critical mass (BLF). Second, galaxies residing in more massive halos today must form, on average, earlier than those inhabiting less massive halos, giving rise to a "downsizing" effect in a ΛCDM universe (e.g. Neistein et al. 2006). Finally, there must be a population of galaxies residing in DM halos with a present-day mass comparable to critical mass that has undergone their formation only recently. Robustly quantifying these expectations is, however, not trivial, as it requires precise knowledge of the assembly history of DM halos and the critical mass for galaxy formation.
In this Letter, we address the last issue. In particular, we apply the recent BLF model (briefly described in Sec. 2.1) together with a high-resolution hydrodynamical cosmological simulation to demonstrate that the existence of a population of late-forming dwarf galaxies that form after redshift z = 3 is a robust cosmological outcome of the ΛCDM model, and highly insensitive to the simulation details, provided stars form in gas that self-gravitates. This prediction stems from the existence of a nontrivial time-dependent critical DM halo mass below which galaxy formation cannot take place and the stochastic growth that characterizes DM halos in ΛCDM, which generally depends on the cosmological parameters (e.g. Lacey & Cole 1993;van den Bosch 2002;Correa et al. 2015, and references therein). We quantify the abundance of late-forming dwarfs expected in a ΛCDM universe, and possible observational counterparts, in Sec. 4.

THE MODEL AND SIMULATION
2.1. The BLF model The BLF model establishes the value of the critical mass for the onset of galaxy formation as a function of time. In this model, galaxy formation takes place in AC halos before CR, and in halos in which gas cannot remain in hydrostatic equilibrium afterward. To determine which halos undergo gravitational collapse to form a galaxy after CR, the model assumes that the gas that falls into dark halos is described by an effective equation of state, which is established by the photoheating background at low densities (n H 10 −4.5 cm −3 ), and by the interplay between photoheating and cooling at high densities. This nontrivial model thus avoids the common assumption that galaxy formation takes place in halos of given (fix) virial 2 temperature, a condition that arises from and Moster et al. (2013), respectively. Galaxies are colored according to their current stellar mass (y-axis).
idealized Jeans arguments. 3 We refer the reader to the original BLF paper for further details and a derivation of the critical mass.

The Simulation
We use a high-resolution hydrodynamical cosmological simulation carried out with the P-Gadget 3 code (last described in Springel 2005) along with the EAGLE model of galaxy formation . The simulation is the same introduced in Benitez-Llambay & Frenk (2020) and we list here only the physical prescriptions relevant for our work. We refer the reader to the original papers for further details. The simulation evolves a periodic cubic volume of side length 20 Mpc, filled with 2×1024 3 DM and gas particles, so that the DM and gas particle mass are m dm = 2.9 × 10 5 M and m gas = 5.4 × 10 4 M , respectively. The adopted Plummer-equivalent gravitational softening is = 195 pc. Gas particles are turned into collisionless star particles once they exceed a density threshold of n thr = 1 cm −3 , a sufficiently high value that ensures that the gas within DM halos becomes self-gravitating before turning into stars. Unlike the original EAGLE simulations, our density threshold for star formation does not depend on metallicity. Our simulation also includes radiative cooling and heating, as tabulated by Wiersma et al. (2009), which in turn includes the Haardt & Madau (2001) UVB radiation field. CR is modeled by turning on the UVB at redshift z re = 11.5 and to ensure that gas is quickly heated to ∼ 10 4 K at z re , an energy of 2 eV per proton mass is instantaneously injected to every gas particle at that time.
DM halos are identified in the simulation using HBT+ (Han et al. 2018), which provides a catalog of "central" and "satellite" DM halos identified in the simulation using a friendsof-friends algorithm with percolation length b = 0.2 in units of the mean interparticle separation, and assigning "bound" particles to each halo based on their binding energies.

Sample Selection
We select simulated galaxies as "central" DM halos that contain more than one stellar particle within their galactic radius, r gal = 0.2 × r 200 , which yields a lower galaxy stellar mass limit of ∼ 5.4 × 10 4 M . As shown in Fig. 1, this criterion imposes a minimum present-day halo mass M 200 ∼ 10 9 M . We also restrict our sample to DM halos with virial mass M 200 10 11 M , as objects above this limit are not of interest for our study because all these systems exceed the AC limit prior to CR (Benitez-Llambay & Frenk 2020), thus making it impossible for these halos to host galaxies that form late. In the selected mass range, the stellar mass of our simulated dwarfs are bounded from above and below by the abundance-matching (AM) constraints from Moster et al. (2013) and Guo et al. (2010), respectively (see Fig. 1). At large halo masses, the mass of the central galaxies are underestimated in our simulation compared to AM expectations, a limitation of the original EAGLE model that does not preclude the analysis that follows.

The Simulated Tail of Late-forming Dwarfs
In order to look for a population of late-forming galaxies, we shall consider the galaxy formation time, t f , defined as the time at which a galaxy first formed a star particle in the simulation. In practice, t f corresponds to the formation time of the oldest star particle found within r gal at z = 0. Fig. 2 shows t f as a function of present-day virial mass for our galaxy sample. Galaxies inhabiting more massive halos at redshift z = 0 form earlier than those hosted by less massive counterparts, although the median t f (shown by the thin red solid line) is weakly dependent on present-day halo mass. More interesting is the fact that the scatter in t f increases significantly at low masses, peaking at about the present-day value of the BLF mass (vertical dashed line). The large scatter at low masses originates from a tail of late-forming dwarfs that we arbitrarily define as galaxies with t f 2.2 Gyr (or z 3), and that constitutes less than ∼ 8% of population of dwarfs with stellar masses M gal 5.4 × 10 5 M in the halo mass range 3 × 10 9 M 200 /M 10 10 . But what determines, in general, the t f versus M 200 relation shown in Fig. 2, and more particularly, what is the origin of the tail of late-forming dwarfs identified in our simulation?
In light of the discussion of Sec. 1, one may speculate that the main trend of the t f versus M 200 relation displayed by our galaxy sample is related to the interplay between the mean assembly history of their host halos and the evolution of the critical mass for galaxy formation. We demonstrate that this  is indeed the case in Fig. 3, which shows t f as a function of halo mass, but as measured at the galaxy formation time rather than at z = 0, M 200,t f . We compare the values measured in the simulation to the time-dependent BLF mass (red dashed line). This critical mass displays a sharp transition from M crit ∼ 3 × 10 7 M to M crit ∼ 10 8 M at z = z re , and results from the fact that galaxy formation is largely limited by atomic hydrogen cooling before CR, and by the ability of the gas to undergo gravitational collapse afterward. Fig. 3 makes it clear that galaxy formation occurs predominantly in halos whose virial mass at t f is comparable to the critical mass. Interestingly, the BLF mass traces the simulated t f versus M 200,t f relation remarkably well, even though neither the BLF model nor the EAGLE model were tuned to do so.

Origin of the Overall Mass Dependence of Galaxy Formation Time
The black dashed lines in Fig. 3 show three mean ΛCDM mass growth histories derived using the Extended Press-Schechter (EPS) formalism (Bond et al. 1991), for halos of present-day mass, 10 11 , 10 10 and 10 9 M . These curves help to clarify the overall mass dependence of the galaxy formation time observed in Fig. 2. Consider, for example, the upper dashed curve, which corresponds to the mean mass growth of a halo with present-day mass, M 200 = 10 9 M . This demonstrates that the bulk population of 10 9 M halos cannot host a galaxy in their center, as the mean mass growth history for these halos is under the BLF mass at all times. On the other hand, massive halos today, those that greatly exceed the present-day critical mass for galaxy formation, have exceeded the BLF mass since before CR. Consequently, these halos host galaxies that have formed very early on, explaining to some extent the observed "downsizing", as argued in Sec. 1. The typical formation time for those intermediatemass halos that eventually exceed the critical mass depends monotonically on halo mass. These arguments indicate that it is possible, in principle, to understand the overall mass-dependent galaxy formation time. To this end, we must compare the mean mass growth of ΛCDM halos to the time-dependent critical mass for galaxy formation, and define t f as the time when the mean halo mass first exceeds the BLF mass. The thick magenta dashed line in Fig. 2 shows the result of this exercise. The agreement between this curve and the median t f that results from the simulation (thin solid line), particularly at high masses, demonstrates that our interpretation is correct. However, this simple exercise results in a divergent t f toward the present-day value of the BLF mass (vertical dashed line), as no galaxy formation can take place below this mass scale, as anticipated. This reasoning would imply that all galaxies inhabiting halos with present-day mass M 200 ∼ M crit,0 are young, contrary to the simulation results, in which only a low fraction of the systems form at late times. It is thus evident that although the mean mass growth of ΛCDM halos is useful to understand the overall trend observed in Fig. 2, it is not sufficient to account for the predominantly old population of galaxies inhabiting halos with present-day mass M 200 M crit,0 . The existence of these galaxies in our simulation is a direct conse-quence of the fact that galaxy formation becomes an increasingly rare event at low halo masses, so the proper modeling of the galaxy formation time at these masses must necessarily take into account the intrinsic scatter in the growth of ΛCDM halos.

Origin of the Scatter in the Galaxy Formation Time
To account for the scatter in the mass growth of halos, we now consider individual (rather than average) EPS growth histories in bins of present-day halo mass in the range, 10 9 M 200 /M 10 11 . We sample the scatter in the assembly of halos by constructing 500 growth histories per mass bin. As in the previous section, we define t f as the time at which the EPS mass first crosses the BLF mass. The thick orange lines of Fig. 2 show the result of this calculation, in particular the median and the 10th and 90th percentiles of the distributions. By considering the intrinsic scatter in the growth history of halos we eliminate the divergence of t f toward the presentday value of the critical mass (vertical line). This is because the great majority of these halos formed their galaxies at earlier times, albeit at much later times than massive halos. This improved model naturally reproduces the increasingly large scatter in t f at low masses. The peak in the 90th percentile around the critical mass is in remarkable agreement with the results from the simulation, so we can safely conclude that the mass-dependent galaxy formation time is simply understood from the interplay between the critical mass for galaxy formation and the time when the halos first exceeded this mass. Within this picture, an unavoidable conclusion is that in ΛCDM there must be a low fraction of halos that have exceeded this critical mass for galaxy formation particularly late (z < 3), thus triggering the first episode of star formation at late times. Late-forming dwarfs can thus be regarded as robust cosmological outcomes, being their formation time largely insensitive to the details of the modeling included in numerical simulations. 4 Finally, DM halos less massive than the critical mass today host predominantly old galaxies, similar to more massive halos. As opposed to late-forming dwarfs, these galaxies populate the low fraction of halos that collapsed unusually early to become more massive than the critical mass at earlier times (see Benitez-Llambay & Frenk 2020, for a detailed discussion on this). Some of these halos would host the observed population of ancient ultrafaint dwarfs (e.g. Simon 2019), which we do not resolve in our simulation.

DISCUSSION AND CONCLUSIONS
Our analysis demonstrates that the existence of a small population of late-forming galaxies, here defined as dwarfs that started forming stars after z < 3, is a robust outcome of the ΛCDM model. This is because their formation depends on the ability of the gas to undergo gravitational collapse, which does not depend on the particularities of the galaxy formation model assumed in simulations. Indeed, Fig. 3 demonstrates that the galaxy formation time is a fair measure of the time when the ΛCDM halos first exceeded the critical mass above which gas cannot remain in hydrostatic equilibrium. Moreover, we showed in Fig. 2 that the BLF model, together with the EPS formalism, enables us to understand not only the dependence of the galaxy formation time on halo mass but also its scatter.
Our results are also robust to the assumed redshift of reionization. Fig. 4 shows the frequency of late-forming dwarfs as measured in our simulation and as obtained from comparing EPS mass growth histories to the BLF mass, assuming that the universe undergoes reionization at z re = 6, 8, 10, and 12. Clearly, the frequency of late-forming dwarfs is largely insensitive to the exact value of z re provided z re 8, a lower limit consistent with recent Planck results (Planck Collaboration et al. 2020).
Finally, we note that our results rely on the idea that galaxy formation largely proceeds in AC halos prior to cosmic reionization (e.g., Bromm & Yoshida 2011, and references therein). Numerical simulations that include physical ingredients missed in our simulation and that are important to address this issue -such as molecular hydrogen cooling and radiative feedback effects-support the idea that the first galaxies indeed form predominantly in atomic cooling halos (e.g. Oh & Haiman 2002;Greif et al. 2008;Wise et al. 2014, and references therein). In particular, the work by Wise et al. (2014) shows that the fraction of halos that host galaxies prior to cosmic reionization vanishes in a narrow range of halo mass below the AC limit. Our results thus rest on assumptions that appear plausible and warrant further scrutiny.

Late-forming Dwarfs in Other Simulations
Interestingly, late-forming dwarfs have already been spotted in cosmological hydrodynamical simulations. For example, using a high-resolution zoom-in simulation of the formation of the Local Group from the CLUES project, Benítez-Llambay et al. (2015) identified two dwarf galaxies with a significant delay in their formation, and they ascribed their origin to the impact of CR. Similarly, using a sample of 15 zoom-in cosmological simulations carried out with the FIRE code, Fitts et al. (2017) identified one dwarf galaxy that formed after z < 1. They, too, ascribed the unusual delayed formation of this galaxy to the effect of CR. Although not discussed by the authors, a small fraction of the isolated low-mass dwarf galaxies analyzed by Garrison-Kimmel et al. (2019) also formed particularly late compared to most dwarfs in their simulations. Finally, recent works have shown that dwarf galaxies that form early can undergo late episodes of star formation once their halos experience significant mergers, a phenomenon related to the fact that DM halos can exceed the critical mass more than once during their lifetime (e.g. Benítez-Llambay et al. 2016;Rey et al. 2020, and references therein). This suggests that the critical mass for galaxy formation does not only establish the onset of galaxy formation for starless halos, but also determines the ability of luminous halos to collect sufficient gas to sustain star formation in their center, and may shape the star formation history of dwarfs. 5 All these results thus strongly support the idea that lateforming dwarfs are rare objects that arise naturally in numerical simulations, regardless of the adopted modeling.

Predicted Number Density
The number density of late-forming dwarfs depends on the abundance of halos that exceed the BLF mass after z 3. In Fig. 5 we show the galaxy stellar mass function of the simulated galaxy population, and of galaxies that formed their first stars after redshift z f . As anticipated throughout our Letter, only a low fraction of galaxies make up the population of late-forming dwarfs (i.e., those with z f < 3). Moreover, the fraction of late-forming dwarfs decreases steadily with decreasing formation redshift and with increasing stellar mass. Finding a massive dwarf (M gal > 10 7 M ) undergoing its formation today becomes thus extremely rare; in fact, dwarfs undergoing their formation at the present day are expected to be faint, with masses comparable to the faint and ultrafaint dwarfs observed nearby (M gal 10 6 M ). Although this particular statement might depend on the particular galaxy formation model included in our simulation, the , make up less than 20% of the simulated galaxy population at the low-mass end, and they contribute even less at larger stellar masses.
good match between the observed galaxy stellar mass function from Baldry et al. (2012) and that measured in our simulation (brown line in Fig. 5) suggests that the stellar masses of the simulated dwarfs are robust, at least in the mass range of overlap.

Observational Counterparts of the Late-forming Dwarfs
It is natural to expect that late-forming galaxies exhibit systematically smaller sizes than the majority of dwarfs. This is because late-forming dwarfs form from gas that dissipates its thermal energy and sinks to the center at late times, preventing their young stars from being stirred up by mergers, as opposed to older dwarfs (see, e.g., Benítez-Llambay et al. 2016, for a discussion of this effect in dwarf galaxies). Our simulation indicates that the half-mass radius of late-forming dwarfs resolved with more than 50 particles is indeed, on average, 30% smaller than that of older dwarfs of similar present-day stellar mass. Besides their compactness, these galaxies form from generally pristine gas. Therefore, the most obvious candidates for the ΛCDM late-forming dwarfs are some of the most metal-poor star-forming blue compact dwarfs (BCDs) known in the Local Volume. BCDs such as I Zwicky 18 or DDO 68 are indeed characterized by compact sizes, ongoing intense star-formation activity in their center, and surprisingly low metallicity. Due to the difficulty of reconstructing star-formation histories with high resolution at early times, it is however unclear whether they contain a substantial population of old stars (e.g. Papaderos et al. 2002;Izotov & Thuan 2004;Pustilnik et al. 2005Pustilnik et al. , 2008Jamet et al. 2010).
Consider, for example, the case of I Zwicky 18. The vigorous star formation at its center at a rate of ∼ 16 × 10 −2 M yr −1 kpc −2 , together with its large gas supplies, low dust abundance (e.g. Wu et al. 2007), low-metallicity of ∼ 1/50 Z (e.g. Aloisi et al. 1999), and slow circular velocity of (38 ± 4.4) km s −1 (Lelli et al. 2012), make this galaxy an ideal analog for the most massive galaxies found in our sample of simulated late-forming dwarfs whose gaseous halo has started the runaway gravitational collapse not so long ago. Also, the extreme properties of I Zwicky 18 have led some authors to argue that this dwarf would be a present-day analog of star-forming galaxies found at high redshift (e.g. Papaderos &Östlin 2012). However, the idea that I Zwicky 18 is a dwarf undergoing its first formation today has been disputed on the grounds that, similarly to other observed BCDs (e.g. Schulte-Ladbeck et al. 2001;Annibali et al. 2003;Vallenari et al. 2005;Aloisi et al. 2005;Makarov et al. 2017), I Zwicky 18 contains RGB stars that put a lower limit constraint to its formation time of about 1 Gyr ago (e.g. Aloisi et al. 2007;Annibali et al. 2013;Sacchi et al. 2016, and references therein). This limit, albeit old in the context of stellar population synthesis models, only constitutes a small fraction of the cosmic scales spanned by the late-forming dwarfs discussed here.
On statistical grounds, however, the scarcity of massive late-forming dwarfs in ΛCDM indicates that most BCDs observed in the local universe cannot be undergoing their formation for the first time today. If our volume is representative of the local universe, our simulation suggests that the great majority of late-forming dwarfs should be much fainter than I Zwicky 18, and much fainter than most BCDs. Indeed, BCDs have stellar masses in the range 10 7 M gal /M 10 9 , whereas Fig. 5 shows that simulated late-forming dwarfs have stellar masses M gal 10 7 M . Whether one of the known BCDs with extreme properties constitutes a true analog of the late-forming dwarfs analyzed here could be answered with deeper photometric observations that reach the oldest main-sequence turnoff (see, e.g., the review by Gallart et al. 2005).
The higher resolution and sensitivity afforded by upcoming observational facilities, in particular the James Webb Space Telescope and the Extremely Large Telescope, will enable precise measurements of resolved star formation histories in these types of galaxies at larger distances, helping to constrain their formation epoch and placing them in the context of the population analyzed in our work. On the other hand, upcoming surveys such as those by the Vera Rubin Observatory may be able to detect ultrafaint dwarfs beyond our Local Group, some of which may well be the late-forming dwarfs discussed here. These endeavors should not be taken lightly, as late-forming dwarfs could provide a powerful avenue to elucidate the past growth of low-mass DM halos, probe a distinctive halo mass-scale, and test the core of our understanding of galaxy formation at the smallest scales.