EDGE: the mass-metallicity relation as a critical test of galaxy formation physics

TL;DR

This paper introduces the EDGE project, studying the formation of small galaxies through high-resolution simulations, revealing how radiative transfer influences star formation regulation and how the stellar mass-metallicity relation constrains galaxy formation physics.

Contribution

It demonstrates the impact of radiative transfer on galaxy self-regulation and highlights the stellar mass-metallicity relation as a key constraint on galaxy formation models.

Findings

Radiative transfer suppresses star formation by keeping gas warm.

Structural galaxy properties are mainly determined by dark matter halos.

Stellar mass-metallicity relation differentiates galaxy formation models.

Abstract

We introduce the "Engineering Dwarfs at Galaxy formation's Edge" (EDGE) project to study the cosmological formation and evolution of the smallest galaxies in the Universe. In this first paper, we explore the effects of resolution and sub-grid physics on a single low mass halo ($M_{\rm halo}=10^{9}~M_\odot$), simulated to redshift $z=0$ at a mass and spatial resolution of $\sim 20~M_\odot$ and $\sim 3$ pc. We consider different star formation prescriptions, supernova feedback strengths and on-the-fly radiative transfer (RT). We show that RT changes the mode of galactic self-regulation at this halo mass, suppressing star formation by causing the interstellar and circumgalactic gas to remain predominantly warm ($\sim 10^4$ K) even before cosmic reionisation. By contrast, without RT, star formation regulation occurs only through starbursts and their associated vigorous galactic outflows. In spite of this difference, the entire simulation suite (with the exception of models without any feedback) matches observed dwarf galaxy sizes, velocity dispersions, $V$-band magnitudes and dynamical mass-to-light-ratios. This is because such structural scaling relations are predominantly set by the host dark matter halo, with the remaining model-to-model variation being smaller than the observational scatter. We find that only the stellar mass-metallicity relation differentiates the galaxy formation models. Explosive feedback ejects more metals from the dwarf, leading to a lower metallicity at a fixed stellar mass. We conclude that the stellar mass-metallicity relation of the very smallest galaxies provides a unique constraint on galaxy formation physics.

Tables

Table 1. Table 1: Photon group energy (frequency) intervals and properties from the ‘Hires+RT’ simulation. The energy intervals defined by the groups are indicated in units of eV by ϵ 0 subscript italic-ϵ 0 \epsilon_{0} and ϵ 1 subscript italic-ϵ 1 \epsilon_{1} . The last four columns show photon properties derived every 10 10 10 coarse time-steps from the stellar luminosity weighted SED model. These properties evolve over time as the stellar populations age, and the mean values are quoted. ϵ ¯ ¯ italic-ϵ \bar{\epsilon} denotes the photon energies, while σ H 2 subscript 𝜎 subscript H 2 \sigma_{\rm H_{2}} , σ H ​ i subscript 𝜎 H i \sigma_{{\rm{H\textsc{i}}}} , σ He ​ i subscript 𝜎 He i \sigma_{{\rm{He\textsc{i}}}} , and σ He ​ ii subscript 𝜎 He ii \sigma_{{\rm{He\textsc{ii}}}} denote the cross sections for ionisation of molecular hydrogen, hydrogen and helium, respectively.

Photon group	${ϵ}_0$ [eV]	${ϵ}_1$ [eV]	$\overline{ϵ}$ [eV]	${\sigma }_{H_2}\left[{\mathrm{cm}}^2\right]$	${\sigma }_{H\text{i}}\left[{\mathrm{cm}}^2\right]$	${\sigma }_{\mathrm{He}\text{i}}\left[{\mathrm{cm}}^2\right]$	${\sigma }_{\mathrm{He}\text{ii}}\left[{\mathrm{cm}}^2\right]$
IR	0.1	1.0	0.6	0	0	0	0
Optical+FUV	1.0	12.0	3.2	0	0	0	0
Lyman-Werner	12.0	13.6	12.6	$1.7\times {10}^{-19}$	0	0	0
UV${}_{H\text{i}}$	13.6	24.59	18.1	$5.1\times {10}^{-18}$	$3.3\times {10}^{-18}$	0	0
UV${}_{\mathrm{He}\text{i}}$	24.59	54.42	35.6	$2.0\times {10}^{-18}$	$5.9\times {10}^{-19}$	$4.3\times {10}^{-18}$	0
UV${}_{\mathrm{He}\text{ii}}$	54.42	$\infty$	64.9	$3.4\times {10}^{-19}$	$8.3\times {10}^{-20}$	$1.2\times {10}^{-18}$	$1.1\times {10}^{-18}$

Table 2. Table 2: Simulations and their z = 0 𝑧 0 z=0 properties. All simulations use the same cosmological zoom initial conditions, are run with RAMSES-RT and reach a minimum cell size of Δ ​ x min = 3 ​ pc Δ subscript 𝑥 min 3 pc \Delta x_{\rm min}=3{\,\rm pc} . The specified quantities are, from left to right: the dark matter particle mass in the deepest refinement region; the gas refinement mass above which new cells will be opened; the total stellar mass formed by z = 0 𝑧 0 z=0 ; the V 𝑉 V -band magnitude computed using sunset ray-tracing; the half-mass radius; the 1D equivalent velocity dispersion; the dynamical mass-to-light ratio within this radius; the enrichment relative to solar; and the mean star formation rate over the time until the galaxy quenches.

Simulation	$m_{\mathrm{DM}}$	$m_{\mathrm{bar}}$	$M_{\star }$	$M_V$	$r_{1/2}$	${\sigma }_{\star }$	$M_{\mathrm{dyn}}/L$	[Fe/H]	${⟨\mathrm{SFR}⟩}^a$
	[$M_{\odot }$]	[$M_{\odot }$]	[${10}^5{M}_{\odot }$]		[pc]	[$\mathrm{km}{s}^{-1}$]	()	[dex]	[$M_{\odot }/$yr]
Fiducial, no feedback	945	161	172.1	-12.0	51	17.6	2	no enrichment	$1.4\times {10}^{-2}$
Fiducial	945	161	1.2	-6.6	313	6.1	287	$-2.65\pm 0.82$	$1.1\times {10}^{-4}$
Fiducial (weak feedback)b	945	161	15.8	-9.4	336	6.1	17	$-1.07\pm 0.68$	$1.4\times {10}^{-3}$
Fiducial, $2\times E_{\mathrm{SN}}$	945	161	1.2	-6.6	280	6.3	212	$-2.66\pm 0.79$	$1.1\times {10}^{-4}$
Fiducial, $10\times E_{\mathrm{SN}}$	945	161	0.52	-5.7	314	6.9	647	$-3.42\pm 0.96$	$4.8\times {10}^{-5}$
Fiducial, $100\times E_{\mathrm{SN}}$	945	161	0.11	-4.0	157	5.6	949	$-4.00\pm 1.20$	$2\times {10}^{-4}$
Fiducial, H2 SF	945	161	1.2	-6.6	285	6.1	189	$-2.49\pm 0.99$	$1.1\times {10}^{-4}$
Fiducial, no H2 physics	945	161	0.53	-5.7	265	5.6	338	$-2.27\pm 1.29$	$4.8\times {10}^{-5}$
Fiducial, HM UV	945	161	0.94	-6.35	283	6.6	287	$-2.59\pm 0.79$	$9.4\times {10}^{-5}$
Fiducial + RT	945	161	0.30	-5.7	308	5.6	398	$-2.31\pm 0.88$	$2.7\times {10}^{-5}$
Fiducial + RT (weak feedback)b	945	161	1.1	-7.1	412	6.5	200.8	$-2.35\pm 0.96$	$9.7\times {10}^{-5}$
Fiducial + RT, H2-based SF	945	161	0.62	-6.5	597	6.5	535	$-1.91\pm 0.65$	$5.3\times {10}^{-5}$
Hires	118	20	2.5	-7.4	370	7.2	171	$-2.61\pm 0.95$	$2.3\times {10}^{-4}$
Hires (weak feedback)b	118	20	20.5	-9.7	311	6.1	12.3	$-1.01\pm 0.91$	$1.8\times {10}^{-3}$
Hires + RT	118	20	0.31	-5.7	203	5.2	226	$-2.6\pm 1.37$	$2.9\times {10}^{-5}$
Hires + RT (weak feedback)b	118	20	2.8	-8.1	368	6.3	67.1	$-1.32\pm 0.85$	$2.3\times {10}^{-4}$
Fiducial, dark matter only	1106	-	-	-	-	-	-
Hires, dark matter only	138	-	-	-	-	-	-

Equations

$${\dot{\rho }}_{\ast }={ϵ}_{\mathrm{ff}}\frac{{\rho }_{\mathrm{g}}}{t_{\mathrm{ff}}}\mathrm{for}{\rho }_{\mathrm{g}}>{\rho }_{\star },$$

$$M_Z=2.09M_{\mathrm{O}}+1.06M_{\mathrm{Fe}}$$

$$M_{\mathrm{dyn}}=3{\sigma }^2{r}_{1/2}/G.$$

$$[Y/X]={\log }_{10}\left(\frac{f_Y/m_Y}{f_X/m_X}\right)-({\log }_{10}{ϵ}_{Y,\odot }-{\log }_{10}{ϵ}_{X,\odot })$$

$$[Fe/H]=\frac{{\sum }_i^N{[Fe/H]}_i{m}_{\star ,i}}{{\sum }_i^N{m}_{\star ,i}}$$

Full text

1 Introduction

In our current $\Lambda$ Cold Dark Matter ($\Lambda$CDM) cosmological paradigm (e.g Komatsu et al., 2011; Planck Collaboration et al., 2014), galaxies form through the successive ‘hierarchical’ mergers of smaller galaxies (e.g. White & Rees, 1978; Fall & Efstathiou, 1980). This theory has been tremendously successful at matching the observed structure in the Universe on large scales (e.g. Springel et al., 2006; Vikhlinin et al., 2009). However, on smaller scales there have been long-standing tensions (e.g. Bullock & Boylan-Kolchin, 2017). In part, such tensions owe to the uncertain mapping between stars and dark matter in low mass dwarf galaxies. This arises because dwarfs are very sensitive to the physics of galaxy formation. Isolated dwarfs can have their star formation suppressed or even extinguished by supernovae (SN) feedback (e.g. Dekel & Silk, 1986; Efstathiou, 2000), reionisation (e.g. Efstathiou, 1992; Barkana & Loeb, 1999; Bullock et al., 2000; Benson et al., 2002), stellar winds and radiative feedback (e.g. Agertz et al., 2013; Hopkins et al., 2014), or some combination of these (e.g. Kravtsov et al., 2004; Read et al., 2006b). If such dwarfs fall into a larger host galaxy, becoming satellites, then they are further buffeted by ram pressure stripping of their interstellar medium (e.g. Grebel et al., 2003; Gatto et al., 2013) and tides (e.g. Kravtsov et al., 2004; Mayer et al., 2006; Read et al., 2006a).

The sensitivity of dwarfs to galaxy formation physics, like the gas density at which stars form (e.g. Kravtsov, 2003; Saitoh et al., 2008) or the details of galactic outflows (e.g. Read et al., 2006b), makes them a natural ‘rosetta stone’ for constraining galaxy formation models. Early work simulating dwarf galaxies focussed on high resolution small box simulations, stopping at high redshift ($z\sim 5-10$) to avoid gravitational collapse on the scale of the box (e.g. Abel et al., 1998; Read et al., 2006b; Gnedin & Kravtsov, 2006; Mashchenko et al., 2008). These simulations demonstrated that stellar winds, supernova feedback and ionising radiation combine to prevent star formation in halo masses below $\sim 10^{7-8}$ M⊙ (Read et al., 2006b; Bland-Hawthorn et al., 2015). Furthermore, once cooling below $10^{4}$ K is permitted and gas is allowed to reach high densities n_{\rm max}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10~{}{\rm cm}^{-3} (requiring a spatial and mass resolution better than \Delta x\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}100 pc and $m_{\rm min}<10^{3}$ M⊙; Pontzen & Governato 2012; Dutton et al. 2018; Benitez-Llambay et al. 2018; Bose et al. 2018), star formation becomes much more stochastic and violent (Mashchenko et al., 2008; Pontzen & Governato, 2012; Dutton et al., 2018). The repeated action of gas cooling and blow-out due to feedback expels dark matter from the galaxy centre, transforming an initially dense dark matter cusp to a core (e.g. Read & Gilmore, 2005; Pontzen & Governato, 2012). Finally, independent of any internal sources of feedback energy, cosmic reionisation can halt star formation in low mass galaxies (e.g. Gnedin & Kravtsov, 2006). Dwarfs that have not reached a mass of M_{200}\lower 2.15277pt\hbox{;\buildrel>\over{\sim};}10^{8} M⊙ (see Section 2.4 for a definition of $M_{200}$) by the redshift that reionisation begins ($z\sim 8-10$; Gnedin & Kaurov 2014; Ocvirk et al. 2018) are gradually starved of fresh cold gas, causing their star formation to shut down by a redshift of $z\sim 4$ (Oñorbe et al., 2015). This is similar to the age of nearby ‘ultra-faint’ dwarf galaxies (UFDs) that have M_{*}\lower 2.15277pt\hbox{;\buildrel<\over{\sim};}10^{5} M⊙, suggesting that at least some of these are likely to be relics from reionisation, inhabiting pre-infall halo masses in the range $M_{\rm 200}\sim 10^{8-9}$ M⊙ (Gnedin & Kravtsov, 2006; Bovill & Ricotti, 2009, 2011; Weisz et al., 2014; Brown et al., 2014; Jethwa et al., 2018; Read & Erkal, 2018).

To compare quantitatively to such very small galaxies in the local Universe (e.g. Simon, 2019), we need simulations with sufficiently high resolution to model the interstellar medium accurately. Most recent work in this area has focused on isolated dwarfs, removing the need to simultaneously capture a large host galaxy like the Milky Way. While much work has recently been done on simulating the smallest dwarf galaxies, the results from different groups, each of whom make different choices for their sub-grid physics model, box size, hydrodynamic and gravity solver, and resolution, are highly discrepant below $M_{\rm 200}\sim 10^{10}$ M⊙ (Oñorbe et al., 2015; Wheeler et al., 2015; Sawala et al., 2015; Read et al., 2016a; Read et al., 2017; Fitts et al., 2017; Munshi et al., 2017; Macciò et al., 2017; Munshi et al., 2018; Revaz & Jablonka, 2018; Wheeler et al., 2018; Smith et al., 2019). While many simulations reproduce the observed structural properties of dwarf galaxies, such as galaxy sizes and velocity dispersions, there are orders of magnitude differences in the stellar mass to halo mass relation, and simulations struggle to reproduce the stellar mass-metallicity relation of the faintest dwarfs (see discussion below; and for a data compilation, see Figures 3 and 6). Furthermore, none of the simulations has yet managed to produce a fully convincing explanation for the existence of star forming ultra faint dwarf galaxies like Leo T (except possibly Wright et al. 2018, see also Verbeke et al. 2015), as cosmic reionisation ought to have evaporated all cold gas from such low mass galaxies. Finally, while observationally there is mounting evidence for dark matter cores in at least some ultra-faint dwarf galaxies (Amorisco, 2017; Contenta et al., 2018; Sanders et al., 2018), some groups find core formation can occur ‘all the way down’ to the very lowest mass dwarfs (Read et al., 2016a; Munshi et al., 2017), while others find that core formation ceases below $\sim 10^{10}$ M⊙ (Oñorbe et al., 2015; Tollet et al., 2016).

In this paper, we introduce a new dwarf galaxy simulation campaign: “Engineering Dwarfs at Galaxy formation’s Edge” (EDGE) with the goal of shedding light on the above discrepancies. There are several new elements to EDGE that make this investigation possible. Firstly, we work at a mass and spatial resolution of $\sim 20\,{\,\rm M}_{\odot}$ and $\sim 3$ pc, respectively, allowing us to better capture the impact of each individual supernova explosion on the forming dwarf. This simplifies the sub-grid modelling of the SN, reducing the need for delayed cooling, momentum capturing schemes, or other similar prescriptions that are required at lower resolution to prevent over-cooling of the SN-heated gas (see e.g. Stinson et al. 2006a; Dalla Vecchia & Schaye 2008; Torrey et al. 2014; Crain et al. 2015 and for a discussion, see Read et al. 2016a). Secondly, we explore the effect of switching on each piece of our sub-grid model for star formation and feedback one at a time, allowing us to assess its role in shaping the final observed properties of the dwarf today. Thirdly, we model isolated dwarfs in a void region down to redshift zero, using a zoom technique. This allows us to explicitly compare our results with observations of ultra-faint dwarfs that can only been seen, at present, in the Local Group. Finally, we set up our cosmological initial conditions using the new genetIC code (Roth et al., 2016; Rey & Pontzen, 2018). This will allow us (in forthcoming papers) to forensically explore the effect of different merger histories and environments on the properties of dwarf galaxies. In this first paper in the series, we present the results of 16 hydrodynamical simulations of a single $M_{\rm 200}=10^{9}{\,\rm M}_{\odot}$ dwarf run at different resolutions and with different sub-grid physics models. In particular, we explore different star formation prescriptions, the effect of SN feedback and the effect of on-the-fly radiative transfer (RT). We emphasise that our goal in running this large suite of simulations is not to ‘calibrate’ our sub-grid physics model. Rather, we seek to address the following questions: (i) which observables are most sensitive to changes in our sub-grid model? and (ii) which physics are most important for regulating star formation and determining the final observed properties of the smallest dwarfs? In further papers in the series, we will explore the role of physics not considered in this work, the role of different merger histories, the effect of increasing halo mass, and the impact of star formation and feedback on the inner density of the dwarf’s dark matter halo.

This paper is organised as follows. In Section 2, we describe the code that we use to run the simulations, the sub-grid physics models that we explore in this work, and how we set up the initial conditions. In Section 3 we present mass growth histories, the stellar mass-halo mass relation and scaling relations for all simulated dwarf galaxies. In Section 4, we compare our findings to observations and previous simulations in the literature and outline the limitations of our models. Finally, in Section 5, we present our conclusions.

2 Simulations

We use Ramses-RT (Rosdahl et al., 2013; Rosdahl & Teyssier, 2015), which is an radiation hydrodynamics (RHD) extension of the cosmological Adaptive Mesh Refinement (AMR) hydrodynamical code Ramses (Teyssier, 2002)111The public code, including all the RHD extensions used here, can be downloaded at: https://bitbucket.org/rteyssie/ramses., to solve the evolution of dark matter (DM), stellar populations, and gas via gravity, hydrodynamics, radiative transfer, and non-equilibrium radiative cooling/heating. For hydrodynamics, we use the HLLC Riemann solver (Toro et al., 1994) and the MinMod slope limiter to construct gas variables at cell interfaces from their cell-centred values. To close the relation between gas pressure and internal energy, we use an ideal gas equation of state with an adiabatic index $\gamma=5/3$. The dynamics of collisionless DM and star particles are evolved with a multi grid particle-mesh solver and cloud-in-cell interpolation (Guillet & Teyssier, 2011). The advection of radiation between cells is solved with a first-order moment method, using the fully local M1 closure for the Eddington tensor (Levermore, 1984) and the Global-Lax-Friedrich flux function for constructing the inter-cell radiation field. With the M1 closure, the collisionless nature of photons is lost and beams are not perfectly maintained (see e.g. Rosdahl et al., 2013). This has the effect that in the case where radiation from many sources is mixed, the radiation flux becomes distorted, by up to a factor two (Decataldo et al., 2019), compared to the real flux which can for example be obtained with ray-tracing methods. These more accurate methods, however, are prohibitively expensive to use in simulations with the number of resolution elements and, more importantly, number of radiation sources that we model in this work. We argue that the modest factor of two errors produced by the M1 closure are likely less significant than several other approximations going into our galaxy formation models.

In the following sub-sections, we describe the set-up of our simulations (multi-frequency radiation, thermochemistry, initial conditions) and the adopted galaxy formation physics (star formation and stellar feedback).

2.1 Radiation

Star particles (see 2.3) are treated as single stellar populations (SSPs) with spectral energy distributions (SEDs) taken from Bruzual & Charlot (2003) (see Rosdahl et al., 2013, for details). We employ six photon groups to account for 1) photoionisation heating, where three groups bracket the ionization energies for Hi, Hei, and Heii, 2) H2 dissociating Lyman-Werner radiation, 3) direct (single scattering) radiation pressure222adopting the “reduced flux approximation” decribed in Appendix B of Rosdahl et al. (2015)., and 4) non-thermal radiation pressure from multi-scattered IR photons. For dust opacities, we adopt $\kappa=10~{}(Z/Z_{\odot})~{}{\rm cm^{2}/g}$ for the IR photons, while for the higher energy photons we assume $\kappa=1000~{}(Z/Z_{\odot})~{}{\rm cm^{2}/g}$, with $Z/Z_{\odot}$ being the gas metallicity in units of the solar value (here taken to be $Z_{\odot}=0.02$). For a full description of our treatment of dust physics, see Rosdahl & Teyssier 2015 and Kimm et al. 2017. Although we include a treatment of radiation pressure in our RT scheme we find that it has negligible impact, likely due to the low dust content in the simulated UFDs; final galaxy stellar masses in two tests without any radiation pressure (not presented here for brevity) end up within $\pm 20\%$ of stellar masses in models including it (but see Wise et al., 2012a).

The group properties (average energies and cross sections to molecular hydrogen, hydrogen and helium) are updated every 10 coarse time-steps from luminosity-weighted averages of the spectra of all stellar populations in the simulation volume, as described in Rosdahl et al. (2013). We do this so that at any time, the cross sections and photon energies are representative of the luminous stellar populations.

In Table 1 we summarize the above information, and present the average group energies and cross sections, over the entire 13.7 Gyr simulation time from one of our high resolution simulations discussed below333Similar values are obtained from all of our simulations where RT is included.. The values only change by a few tens of percent over the course of a simulation as the contributions from both old and metal-rich stellar populations increase.

2.2 Gas thermochemistry

For $T>10^{4}$ K, the contribution from metals to gas cooling is computed using tables generated with Cloudy (Ferland et al., 1998, version 6.02). We adopt the UV background of Faucher-Giguère et al. (2009) as our standard setting, but also consider that of Haardt & Madau (1996) (see Section  2.5). The homogenous UV background is necessary to include as we are not capturing radiation sources outside the zoom region. It is not treated using radiative transfer but in the commonly adopted optically thin, cell-by-cell heating approximation (e.g. Katz et al., 1996). Following Rosdahl & Blaizot (2012) self-shielding against the homogenous UV background is modelled by applying a damping function to the photoionisation rate, $\Gamma_{\rm ss}=\Gamma_{\rm UV}\exp(-n_{\rm H}/10^{-2}~{}{\rm cm^{-3}})$, for hydrogen densities $n_{\rm H}>10^{-2}~{}{\rm cm^{-3}}$. Self-shielding against radiation modeled using RT is treated self-consistently.

For $T\leq 10^{4}$ K, we use the fine structure cooling rates from Rosen & Bregman (1995). The non-equilibrium hydrogen and helium thermochemistry, coupled with the local ionising radiation, is performed with the quasi-implicit method described in Rosdahl et al. (2013) via photoionisation, collisional ionisation, collisional excitation, recombination, bremsstrahlung, homogeneous Compton cooling/heating off the cosmic-microwave background, and di-electronic recombination. We account for the formation, advection, destruction and cooling of molecular hydrogen (H2) (see Nickerson et al., 2018, for details), and its coupling to the radiation field, in all simulations unless otherwise stated.

Along with the temperature, and photon fluxes, we track, in every cell, the fractions of neutral hydrogen, ionised hydrogen, singly, and doubly ionised helium ($x_{\rm H\textsc{i}}$, $x_{\rm H\textsc{ii}}$, $x_{\rm{He\textsc{ii}}}$, $x_{\rm{He\textsc{iii}}}$, respectively), and advect them with the gas as passive scalars. The thermochemistry is operator split from the advection of gas and radiation and performed with adaptive time step sub-cycling on every RT time-step. To keep the computational costs low, we use a reduced speed of light, $\bar{c}=10^{-2}c$ (Rosdahl et al., 2013), where c is the true speed of light.

2.3 Galaxy formation physics

Star formation follows a Schmidt law,

[TABLE]

where $\rho_{\rm g}$ is the gas density, $\rho_{\star}$ the density threshold of star formation, $t_{\rm ff}=\sqrt{3\pi/32G\rho}$ is the local gas free-fall time and $\epsilon_{\rm ff}$ is the star formation efficiency per free-fall time. We adopt $\rho_{\star}=300~{}{m_{\rm H}}~{}{\,\rm{cm^{-3}}}$, and sample Eq. 1 stochastically on a cell-by-cell basis at every fine simulation time step using $300{\,\rm M}_{\odot}$ star particles444Chosen to accommodate the (discrete and stochastic) mass loss from several SN type II explosions, which is not the case if initial star particle masses are set too close to the mass of individual massive stars in our current feedback scheme. Note also that star particle masses,on average, are reduced by up to $50\%$ due to stellar evolution (e.g. Leitner, 2012). (see Agertz et al., 2013). Furthermore, we only allow stars to form from cold gas ($T<100$ K). In a subset of simulations, we require star formation to only occur in molecular gas, as motivated by the observed close to linear relation between $\Sigma_{\rm mol}$ and $\Sigma_{\rm SFR}$ (e.g. Bigiel et al., 2008). In those simulations, $\rho_{\rm g}\rightarrow f_{{\rm H}_{2}}\rho_{\rm g}$ in Eq. 1, where $f_{{\rm H}_{2}}$ is the molecular hydrogen fraction in a cell (see also Gnedin et al., 2009; Gnedin & Kravtsov, 2010; Christensen et al., 2014). We note that while a correlation between molecular hydrogen and star formation is well motivated both theoretically and empirically, it is not well established in the low metallicity regimes probed in this work555The cooling and star formation time scales can in metal poor environments be shorter than the time scale for reaching an equilibrium chemical state for which the gas would be H2-dominated. Star formation can under such conditions correlate with atomic gas (e.g. Krumholz, 2012). The non-equilibrium treatment of H2 in our adopted method (Nickerson et al., 2018) mitigates this issue (see also Krumholz & Gnedin, 2011). ([Fe/H]$\lesssim-2$, Glover & Clark, 2012; Krumholz, 2012). For this reason we consider the use of H2-based star formation models in the UFD regime as exploratory (see also Munshi et al., 2018).

Observationally, $\epsilon_{\rm ff}$ averages $1\%$ on galactic kpc scales (Bigiel et al., 2008) as well as in Milky Way giant molecular clouds (GMCs) (Krumholz & Tan, 2007), albeit with a spread of several dex (Murray, 2011; Lee et al., 2016). Recently Grisdale et al. (2019) demonstrated how high efficiencies ($\epsilon_{\rm ff}\sim 10\%$) on scales of parsecs, coupled to a feedback budget like the one adopted here, provides a close match to the observed (i.e. emerging) efficiencies on scales of individual GMCs. Motivated by these findings we adopt $\epsilon_{\rm ff}=10\%$.

We adopt the stellar feedback budget described in Agertz et al. (2013). Briefly, this feedback prescription includes the injection of energy, momentum, mass and heavy elements over time from SNII and SNIa explosions and stellar winds into the surrounding ISM. In contrast to Agertz & Kravtsov (2015), we do not adopt a subgrid model for radiation pressure, as this is self-consistently treated by the RT solver (see Section  2.1, and Section  2.5 for which simulations adopt RT). Each mechanism depends on the stellar age, mass and gas/stellar metallicity, calibrated on the stellar evolution code STARBURST99 (Leitherer et al., 1999), treating each formed stellar particle as a single-age stellar population with a Kroupa (2001) initial mass function (IMF). Feedback is done continuously at the appropriate times when each feedback process is known to operate, taking into account the lifetime of stars of different masses within a stellar population.

We track iron (Fe) and oxygen (O) abundances separately, and advect them as passive scalars. When computing the gas cooling rate, which is a function of total metallicity, we construct a total metal mass as:

[TABLE]

according to the solar abundances of alpha (C, N, O, Ne, Mg, Si, S) and iron (Fe, Ni) group elements of Asplund et al. (2009).

SNe explosions are modelled as discrete events, and we follow the approach by Kim & Ostriker (2015) (see also Martizzi et al., 2015) and inject the full momentum generated during the Sedov-Taylor phase if a supernova cooling radius is not captured with at least 6 grid cells, otherwise we inject a $E_{\rm SN}=10^{51}$ ergs of thermal energy (see Agertz et al., 2015, for details) and allow for the hydrodynamic solver to track the buildup of momentum. At the numerical resolution adopted here, $>90\%$ of all SN explosions are resolved by at least 6 cooling radii in our fiducial simulation. We note that we do not enforce any additional refinement criterion to achieve this; the (Lagrangian) mass based refinement scheme, discussed in Section 2.5, is enough for this to be satisfied. This feedback budget has been shown to lead to Milky Way disc galaxies in close agreement to observations (Agertz & Kravtsov, 2016, Agertz et al. in prep), bursty star formation and realistic properties of dwarf galaxies (Read et al., 2016a; Read et al., 2016b), and an interstellar medium (ISM) and structure of giant molecular clouds in excellent agreement with observations (Grisdale et al., 2017, 2018).

Finally, in order to account for enrichment from unresolved Population III (Pop III) star formation, we adopt a pre-existing metal floor at $Z=10^{-3}Z_{\odot}$ (e.g. Wise et al., 2012b; Jaacks et al., 2018) added to the oxygen field666Motivated by the abundance ratios for alpha group elements observed in extremely metal poor stars (Iwamoto et al., 2005), but we note that Pop III yields are extremely uncertain, and that abundance ratios and the overall metal content is likely dependent on environment (e.g. Jaacks et al., 2018).. We note however that H2 cooling is the main coolant at $Z\lesssim 10^{-2}Z_{\odot}$, which is the regime we are exploring in this work. As soon as star formation and enrichment begins, metal line cooling, and subsequently atomic line cooling, also become important coolants (Wise et al., 2014). The inclusion of a Pop III floor therefore has a small effect on star formation properties in our simulations, and tests with $Z\leq 10^{-4}Z_{\odot}$ show that our choices have no impact on any of the conclusions presented in this paper.

2.4 Initial conditions

In this work we study the cosmological formation of dwarf galaxies forming in isolation, hence removing complexities such as gas stripping during infall and environmental star formation quenching. Being isolated, deformation due to tides is also minimized, and we can study the “pristine” galaxy formation scenario in a $10^{9}{\,\rm M}_{\odot}$ dark matter halo. To generate initial conditions, we use the code genetIC (Roth et al., 2016; Rey & Pontzen, 2018) which will in future work allow us to explore a continuum of alternative merger histories for the same galaxy. For the present paper, the most important capability of genetIC is simply to recursively refine regions of the simulation on extremely fine grids. We first generated initial conditions for a simulation with box size of $L_{\rm box}=50~{}{\,\rm Mpc}$ with cosmological parameters $\Omega_{\rm m}=0.309$, $\Omega_{\Lambda}=0.691$, $\Omega_{\rm b}=0.045$ and $H_{0}=67.77~{}\rm{km}\,\rm{s}^{-1}~{}{\,\rm Mpc}^{-1}$, in line with data from the PLANCK satellite (Planck Collaboration et al., 2014).

We simulated this volume with only dark matter from $z=99$ to $z=0$, using $512^{3}$ resolution elements (giving a particle mass of $m_{\rm DM}=3.8\times 10^{7}{\,\rm M}_{\odot}$). Next, we picked the largest void and resimulated it (again with only dark matter) at the equivalent of $2048^{3}$ resolution ($m_{\rm DM}=4.9\times 10^{5}{\,\rm M}_{\odot}$), adding the appropriate small scale power to this grid. We then identified dark matter haloes within the void at $z=0$ using the HOP halo finder (Eisenstein & Hut, 1998) and computed their virial masses, $M_{200}$, defined as the mass inside of a spherical volume encompassing $200$ times the cosmic critical density $\rho_{\rm crit}=3H(z)^{2}/8\pi G$, where $H(z)$ is the Hubble parameter. The radial extent of this spherical volume defines the virial radius, $r_{200}$.

To find haloes in isolation, we computed pairwise distances, measured from the edges of their individual virial radii. For halo $n$ and $m$, the distance is $D_{n,m}=|\mathbf{r}_{n}-\mathbf{r}_{m}|-r_{200,n}-r_{200,m}$, where $\mathbf{r}$ is the position of a halo center. Expressed in units of the virial radius of halo $n$, $I_{n,m}=D_{n,m}/r_{200,n}$. As a measure of isolation for halo $n$, we define an isolation parameter as the minimum of $I_{n,m}$, i.e. $I_{n}={\rm min}(I_{n,m})$. At $z=0$, we only considered halos in the mass range $0.8\times 10^{9}{\,\rm M}_{\odot}<M_{200}<1.2\times 10^{9}{\,\rm M}_{\odot}$, with $I_{n}>10$ for halo pairs with $M_{200,m}/M_{200,n}>1$.

From the filtered halo catalogue, we picked a halo that at $z=0$ had a mass of $M_{200}\sim 10^{9}{\,\rm M}_{\odot}$ and was visually isolated from massive dark matter filaments. Having found a dark matter halo satisfying our selection criteria, we identified all particles belonging to this halo out to $2\times r_{200}$ at $z=0$ and traced the particles back to the initial conditions ($z=99$). For the Lagrangian region of this halo, we generated separate new initial conditions at the equivalent of $16\,384^{3}$ resolution ($m_{\rm DM}=1106~{}{\,\rm M}_{\odot}$), again adding small scale density fluctuations compatible with the background. These ‘fiducial’ initial conditions are then modified to include baryons as well as dark matter, forming the basis for the simulation suite described below. For the ‘Hires’ simulations, we further refined the Lagrangian region by a factor of 8 in mass.

2.5 Simulation suite

For halos of a given mass, there will be a diversity of cosmologically-determined mass accretion histories and environments giving rise to a spread in final observed properties at $z=0$. Our EDGE project will ultimately explore this diversity using “genetic modification” (Roth et al., 2016; Rey & Pontzen, 2018). However, before assessing the connection between history and observables, we need a firm handle on theoretical uncertainties resulting from the small-scale star formation and feedback physics. While there will never be a ‘complete’ account of feedback physics, understanding the leading-order uncertainties and their implications for interpreting future simulations is a key first step.

This first work therefore probes the effect of differing physics implementations, and we undertake this study at two numerical resolutions. In the high resolution simulations, the dark matter particle resolution is $m_{\rm DM}=118{\,\rm M}_{\odot}$, and the equivalent baryon resolution (the mass of baryons at the finest grid level in the initial conditions) is $m_{\rm bar}=\Omega_{\rm b}/\Omega_{\rm m}m_{\rm DM}=20{\,\rm M}_{\odot}$. For our fiducial resolution, the corresponding numbers are $m_{\rm DM}=945{\,\rm M}_{\odot}$ and $m_{\rm bar}=161{\,\rm M}_{\odot}$. In all simulations we reach a mean physical resolution of $\Delta x=3~{}{\,\rm pc}$ at all times in the inner parts of dark matter haloes and in the ISM. Refinement is based on a pseudo-Lagrangian approach, where a cell is split if its mass $m_{\rm cell}$ exceeds $8\times m_{\rm bar}$, where the cell mass accounts for both stars and gas. In addition, a cell is allowed to refine if it contains more than 8 dark matter particles. All simulations are run from $a=0.01$ ($z=99$) to $a=1$ ($z=0$), and simulation snapshots are stored every $\Delta a=0.01$.

For both resolution settings, we adopt the standard galaxy formation physics presented above and refer to these as ‘Fiducial+RT’ and ‘Hires+RT’, and ‘Fiducial’ and ‘Hires’ when RT is not included. We reemphasise that without on-the-fly RT, we do not consider any subgrid model of radiative feedback (e.g. Agertz et al., 2013), only feedback from stellar winds and supernovae. For the fiducial resolution setting, we carry out a suite of simulations where the sensitivity to galaxy formation physics is tested as follows: increasing the strength of SN feedback ($E_{\rm SN}=2,10$ and $100\times 10^{51}$), not considering H2 physics and the associated gas cooling, changing the UV background field (changing from the fiducial Faucher-Giguère et al. 2009 UV background to that of Haardt & Madau 1996), adopting an H2 based star formation model (see Section  2.3) and studying the impact of coupling it to RT. In the latter model, molecular hydrogen can be dissociated by Lyman-Werner radiation, forcing stars to only form in environments where H2 is self-shielded, as well as shielded by dust.

For both resolution settings, with and without RT, we also carry out simulations where we artificially limit the efficiency of SN wind driving by limiting the maximum allowed temperatures of the gas to $T_{\rm max}=10^{8}$ K and restrict the maximum allowed velocities of feedback, upon injection, to \mbox{v}_{\rm fb,max}=1000~{}\rm{km}\,\rm{s}^{-1} (approximately leading to post-shock temperatures of $\sim 10^{8}$ K)777In the adopted feedback model, supernova ejecta from $8~{}{\,\rm M}_{\odot}$ stars travel with velocities of $\sim 3000~{}\rm{km}\,\rm{s}^{-1}$. We refer to these models as ‘weak feedback’. Although the exact values of such feedback limiters are arbitrary, these tests illustrate the role of hot, fast-moving winds in regulating ultra-faint dwarf formation, as we will demonstrate below. We summarise the entire simulations suite, as well as the $z=0$ properties of the central dwarf galaxy, in Table 2.

3 Results

Figure 1 shows the gas surface density, temperature and metallicity in a $5\times 5{\,\rm kpc}^{2}$ region centered on the main dwarf galaxy in the ‘Hires’ and ‘Hires+RT’ simulations at $z=10,8,6.5$ and $5$. The maps illustrate the dramatic effect of stellar feedback, and how the inclusion of radiative feedback changes the mode of galaxy formation. In ‘Hires’, star formation leads to large-scale hot ($T>10^{6}\,{\rm K}$) SN driven outflows at all times before $z\sim 5$. The enriched winds efficiently mix into the intergalactic medium (IGM), leading to an average background [Fe/H]$>-3$ even at early times ($z\sim 10$).

This is in stark contrast to the RT counterpart, where the gas is kept warm ($T\sim 10^{4}$ K) and large scale metal rich outflows are relatively weak until $z\lesssim 6.5$. Before this epoch, the IGM is metal poor ([Fe/H]$\lesssim-4$), and metal line cooling is subdominant to H2 cooling. At $z=5$, reionisation has been operating for $\sim 500$ Myr, and large scale filamentary structures have been evaporated, with stars still forming from the residual cold gas from earlier accretion epochs. At this time, both simulations feature a circumgalactic medium enriched to [Fe/H]$\sim-2$. ‘Hires’, due to its early intense outflows, retains an enriched intergalactic medium to larger radii than ‘Hires + RT’, and the RT simulation also features a visibly denser circum-galactic medium compared to ’Hires’. In the next sections we identify differences in the interstellar medium that give rise to the modified outflow behaviours.

3.1 Mass growth histories

Figure 2 shows the buildup of stellar and dark matter halo masses ($M_{200}$) for the simulations in Table 2. Stellar masses ($M_{\star}$) are throughout this paper defined as the total mass of stars in the inner part of the dark matter halo ($r<0.25R_{200}$). The left hand panel includes the models targeting radiative transfer and H2 physics, while the right panel focuses on variations in feedback strength and numerical resolution. Star formation starts at $z\sim 12-13$ in all models, with delays introduced by variations in how H2 physics is treated; see below. The virial temperature of the halo at $z<10$ is $T_{\rm vir}\sim 5000-9000$ K, i.e. below the ‘atomic-cooling’ regime ($T_{\rm vir}\gtrsim 10^{4}$ K). As reionisation heats the gas to $T\sim 10^{4}$ K, this leads to a cessation of gas accretion and ultimately quenching of star formation. Star formation from the already cold and dense (self-shielded) gas in the ISM can, however, proceed throughout the epoch of reionisation. Indeed, the last stars form around $z\sim 4-5$ in all simulations.

Without any stellar feedback, the final stellar mass is close to $M_{\star}=2\times 10^{7}{\,\rm M}_{\odot}$. By including stellar feedback (but neglecting RT), stellar masses are lowered by over two orders of magnitude to $M_{\star}=1.2\times 10^{5}{\,\rm M}_{\odot}$ in ‘Fiducial’, and $\sim 2.5\times 10^{5}{\,\rm M}_{\odot}$ in ‘Hires’. This result is independent of the adopted UV background, with ‘Fiducial, HM UV’ and ‘Fiducial’ having very similar mass growth histories. We next turn to how adopted galaxy formation physics affects stellar masses using our fiducial simulations suite.

3.1.1 Impact of molecular hydrogen

All of our simulations include H2 formation, destruction and cooling (see Section 2.2), with the exception of the test run ‘Fiducial, no H2 physics’. However typically we determine star formation based on local density and temperature of all gas (Section 2.3). In one variant (‘Fiducial, H2 SF’) we tie the star formation exclusively to the molecular component. These two variants thus allow us to determine the overall importance of molecules in our formulation.

In order to reach sufficiently high H2 gas fractions for star formation, H2 needs to be self-shielded. In the low metallicity ISM of the simulated ultra-faint dwarfs (see Section 3.4), this requires high densities, and we find that star formation occurs only at densities $n\sim 10^{3}-10^{4}~{}{\rm cm}^{-3}$ in ‘Fiducial, H2 SF’. As a result, the onset of star formation is slightly delayed ($z\sim 10$) compared to the fiducial model; nonetheless, as soon as star formation starts, the stellar mass growth quickly catches up with the ‘Fiducial’ galaxy. Overall, therefore, the star formation law is not strongly sensitive to molecules in itself. This insensitivity may not comes as a surprise as star formation anyway is restricted to dense ($n>300~{}{\rm cm}^{-3}$) and cold gas ($T<100$ K), which benefits molecular hydrogen formation (see also Hopkins et al., 2018).

On the other hand, completely neglecting H2 cooling (‘Fiducial, no H2 physics’) delays the onset of star formation more significantly (to $z\sim 8$); the final stellar mass is then suppressed by a factor of two relative to ‘Fiducial’. This indicates that the onset of star formation in cold dense gas is sensitive to details of the cooling function. Once star formation starts, enrichment and metal line cooling allows approximately the same gas reservoir to be available for star formation in all models after $z\simeq 8$, limiting the impact on the final stellar content. We conclude that molecular hydrogen physics is a subdominant uncertainty relative to other factors that we explore below.

3.1.2 Impact of supernova feedback strength

Increasing the available energy per SN explosion by a factor of two (‘Fiducial, $2\times E_{\rm SN}$’) has almost no effect on the stellar mass growth histories. However, an increase by a factor 10 (or 100) results in a suppression by a factor of 2 (or 10) in final stellar masses. This reveals a weak, sub-linear dependence of the amount of star formation on the strength of SN feedback, compatible with results for simulated disc galaxies (Benincasa et al., 2016).

Limiting the effect of stellar feedback (‘Fiducial, weak feedback’), by artificially imposing numerical ceilings on the allowed supernova gas temperature $T_{\mathrm{max}}=10^{8}\,\mathrm{K}$ and velocities ($v_{\mathrm{fb,max}}=10^{3}\,\mathrm{km\,s^{-1}}$), results in an increased stellar mass by an order of magnitude. This stems from the fact that a large fraction of supernova explosions occur in low density gas where SN bubble temperatures and velocities exceed the ceilings (see Section 2.5). Capturing this extreme temperature gas phase, and the associated fast ejecta velocities, is thus essential to capture the full effects of energetic feedback at high resolution.

3.1.3 Impact of radiative feedback

As previously discussed with reference to Figure 1, including radiative feedback (‘RT’ simulations) generates a major shift in the early behaviour of our galaxy. The ‘Fiducial’ and ‘Hires’ models with RT all result in a final galaxy with stellar masses $M_{\star}\approx 3\times 10^{4}{\,\rm M}_{\odot}$ – an additional factor of $\sim 5-10$ reduction from pure SN regulation, and the lowest stellar masses recovered in our simulation suite. Radiative feedback coupled to H2 based star formation, i.e. allowing for H2 destruction by Lyman-Werner radiation from young stars (‘Fiducial + RT, H2-based SF’), particularly suppresses star formation rates at early times ($z\gtrsim 6$), although final stellar masses are close to that of ‘Fiducial + RT’.

Without feedback, star formation rates (averaged in age bins of width $100$ Myr) reach $\dot{M_{\star}}\sim 10^{-2}\rm M_{\odot}\,yr^{-1}$ at a lookback time of $t_{\rm lookback}\sim 12.5-13\,{\rm Gyr}$. The order of magnitude change in final stellar mass when feedback processes are introduced are mirrored in the average star formation rates, with ‘Fiducial’ and ‘Fiducial+RT’ models peaking at $\sim 10^{-4}\rm M_{\odot}\,yr^{-1}$ and $\sim{\rm few}\times 10^{-5}\rm M_{\odot}\,yr^{-1}$ respectively, when averaged over $100$ Myr windows.

It is worth noting that this effect is far greater than the differences generated by changing numerical resolution. ‘Hires + RT’ forms a total of $3.1\times 10^{4}\,{\,\rm M}_{\odot}$ compared to $3.0\times 10^{4}$ for ‘Fiducial + RT’. Enabling RT actually seems to minimise resolution sensitivity, most likely because more gas is kept in a warm and relatively diffuse phase.

3.2 The stellar mass – halo mass relation

In Figure 3, we show the relation between the galaxy stellar mass and dark matter halo mass (the $M_{\star}-M_{200}$ relation) for our simulations. For clarity, we show results from the ‘Hires’ models; the corresponding ‘Fiducial’ runs show the same trends, and Table 2 gives data on the entire simulation suite. The left panel shows a comparison to observations, whereas the right panel focuses on comparisons with existing simulation suites. In both cases we also show the results from abundance matching (AM) as implemented by Read et al. (2017), who used the ‘field’ galaxy stellar mass function from the Sloan Digital Sky Survey (SDSS) and the halo mass function from the cold dark matter Bolshoi simulation (see also Behroozi et al., 2013). Below $M_{200}\sim 5\times 10^{9}{\,\rm M}_{\odot}$ ($M_{\star}\sim 10^{7}{\,\rm M}_{\odot}$), the relation assumes a power-law extrapolation of the SDSS stellar mass function, as indicated by the dashed lines.

Note that the simulated stellar masses are taken ‘as is’, i.e. we do not account for any uncertainties from observational colour fitting procedures. Such uncertainties could lead to observational underestimates of a factor up to $\sim 2$ (e.g. Munshi et al., 2013). Moreover, the simulated $M_{200}$ is the mass of dark matter - we do not include the baryons which at $z=0$ are negligible.

The left panel data compilation consists of local group dwarf galaxies with estimated dynamical masses (Read et al., 2017; Read & Erkal, 2018), including galaxies thought to form in halo masses compatible with our simulated galaxies ($M_{200}\sim 10^{9}{\,\rm M}_{\odot}$); LeoT, Eridanus-II and Carina. As discussed by Read & Erkal (2018), the gas rich, star-forming dwarf irregulars have systematically higher estimated galaxy formation efficiencies ($M_{\star}/M_{200}$) compared to the (non star-forming) satellites. How some low-mass galaxies manage to be star-forming at $z=0$ is not yet clear; see Wright et al. (2018) for a possible line of explanation. In any case, given that our simulated galaxies are quenched reionisation fossils, a direct comparison to the quenched satellites is more appropriate.

The form of the $M_{200}-M_{\star}$ relation in this low mass regime is highly uncertain, but one can immediately rule out the ‘No feedback’ case as overforming stars by at least one dex. We find that the suppression of stellar mass from supernova feedback (‘Hires’, and also in the ‘Fiducial’ case which is not plotted) brings galaxies close to AM predictions and local group dwarfs, with $M_{\star}/M_{200}\sim 2.5\times 10^{-4}$). Almost an additional order of magnitude suppression comes from radiative feedback (‘Hires + RT’), leading to $M_{\star}/M_{200}\sim 4\times 10^{-5}$), below the AM extrapolation, but within the uncertainties on Eridanus-II which is the smallest UFD in our observational sample. The model with weak stellar feedback predicts an efficiency of $M_{\star}/M_{200}\sim 2\times 10^{-3}$ which is in broad agreement with AM and dIrrs — but in tension with quenched UFDs, which are the most appropriate point of comparison.

More observational data on this relationship and its scatter would be helpful to make clearer comparisons to simulations. From the theoretical side, a comparison to existing simulations in the literature is shown in the right panel of Figure 3, highlighting that different studies currently make very different predictions for the light-to-mass ratio in the faintest objects. We will discuss this point further in Section 4.1.2, but it should already be clear from the Figure that there is little agreement on the overall effects of feedback in this regime.

In summary, star formation rates in low mass dwarf galaxies are extremely sensitive to the detailed feedback physics, with radiative feedback contributing significantly to suppression of star formation. Next we demonstrate how these differences arise by studying the density and thermal structure of the ISM.

3.3 Structural differences of the interstellar medium

In Figure 4, we show, for a subset of our simulations in the fiducial suite, the mass weighted density (left) and the temperature probability distribution functions (PDFs; right) of the gas in the inner parts of the dark matter halo ($r<0.25R_{200}$). As the ISM is highly dynamical in low mass galaxies, a single snapshot in time cannot represent the average state of the galaxies. To alleviate this issue888We note that a more robust analysis, not affected by limited time resolution, would have required mass-weighted PDFs to be computed on-the-fly at every simulation time step., the PDFs are mass-weighted averages of all simulation outputs in the redshift range $4<z<13$.

Without feedback (grey PDFs), the majority of the ISM’s mass is locked up in a population of dense and cold ($T\approx$10 K) star forming clouds reaching densities as high as $n\sim 10^{4}~{}{\rm cm}^{-3}$. Enabling SN feedback (black PDF) suppresses the fraction of mass in this phase, and reduces the maximum densities by a factor of 100. Note that most gas, by mass, is still retained in the cold gas phase at $T<100~{}$K (below our adopted star formation temperature threshold).

Adding radiative feedback (blue PDFs) has a dramatic effect on the temperature structure, with the peak of the temperature PDF shifted from $T\approx 10~{}$K to $T\approx 10^{4}~{}$K. The large effect of RT on the thermodynamical state of the ISM, and resulting change in the gas mass available for star formation is the primary mechanism for suppressing galactic star formation rates in these models. The suppression is achieved in a relatively ‘gentle’ fashion, without needing to blow the entire interstellar medium from the galaxy. In fact, we find that before reionisation the total gas mass within the galaxy in ‘Fiducial + RT’ is almost identical to the ‘No feedback’ case, despite the two simulations being at opposite extremes in terms of their star formation rates. The total mass in baryons (gas and stars) in the inner halo ($r<0.25R_{200}$) differ however, with ‘No feedback’ having over 3 times as much baryonic mass at $z\sim 6$ compared to ‘Fiducial + RT’, and almost 10 times the baryonic mass of ‘Fiducial’.

The gentle suppression of star formation in the ‘Fiducial + RT’ simulation is also evident from the gas density PDF, as it features more dense gas than the other models, except for ‘Fiducial, no feedback’. Thus star formation is being suppressed without destroying dense clumps. As discussed in Section 1, radiative feedback operates immediately upon the birth of massive stars, which enables the remaining dense cloud to be gently heated. By way of contrast, the first SN in a stellar population only explodes after $\sim 4$ Myr, which corresponds to several free-fall times in the dense star forming gas999Note that stellar winds operate before the first SN explosions in all models, but the low metallicity makes them inefficient at regulating star formation.. These very different regulation modes will almost certainly have implications for the predicted dark matter distributions in UFDs which we will study in future work; see also Section 4.

Having established the sensitivity of different sub-grid models and numerical resolution on the star formation rates for this isolated dwarf, in the next section we confront our suite of simulations with observations of nearby dwarf galaxies at $z=0$.

3.4 Dwarf galaxy scaling relations

In Figure 5, we compare simulated and observed dwarf galaxy properties. We focus on $z=0$ relations between $V$-band magnitudes ($M_{V}$), half-mass radii ($r_{1/2}$), stellar velocity dispersions ($\sigma_{\star}$) and dynamical mass-to-light ratios ($M_{\rm dyn}/L$). The results use the same point styles for the four simulations already shown in Figure 3; the remaining simulations are shown as small circles. Simulations with RT switched on are highlighted in blue; the ‘Fiducial weak feedback’ simulation is highlighted in yellow; the ‘Fiducial $10\times E_{\mathrm{SN}}$’ is shown in red; the ‘Fiducial’ simulation is shown in black; and all other simulations in the suite are shown in grey. All results are also reported in Table 2. Observational data are compiled from McConnachie (2012), Kirby et al. (2013), Kirby et al. (2014), and the recent compilation for UFDs by Simon (2019)101010These references contain overlapping sources, sometimes with slightly discrepant values for the quantities under consideration. When data is found to overlap, the latest compilation in Simon (2019) take precedence, and for additional sources the data in Kirby et al. (2013) and Kirby et al. (2014) supersede the data in McConnachie (2012).

$V$-band magnitudes are computed using all stars in the main dwarf galaxy halo using the SUNSET111111Publicly available as a part of the RAMSES distribution. ray-tracing code, assuming a Kroupa IMF (Kroupa, 2001). Velocity dispersions are 1D equivalents, i.e. $\sigma_{\star}=\sigma_{\rm\star,3D}/\sqrt{3}$, where $\sigma_{\rm\star,3D}=\sqrt{\sigma_{\rm\star,x}^{2}+\sigma_{\rm\star,y}^{2}+\sigma_{\rm\star,z}^{2}}$. To allow for closer comparison of the dynamical mass-to-light ratios to the observational data, we follow Kirby et al. (2014) and estimate the dynamical mass using the simple mass estimator:

[TABLE]

As can be seen in Figure 5, all simulations, with exception of the outlier ‘Fiducial, no feedback’, have global galaxy properties compatible with observations. In this set of simulations, $M_{V}$ lies predominantly in the range -5.5 to -6.5, with the simulations including RT being the least luminous. Models with ‘weak feedback’ reach $M_{V}\sim-9.5$. Despite the range in $M_{V}$, $z=0$ galaxy sizes and velocity dispersions are found to be quite insensitive to the different galaxy formation scenarios, with $r_{1/2}\sim$ a few 100 pc and $\sigma_{\star}\sim 5.5-6.5~{}\rm{km}\,\rm{s}^{-1}$ for all galaxies. This weak dependence on final galaxy masses reflects the dominant contribution of the dark matter halo to the gravitational potential. Indeed, dynamical mass-to-light ratios for all galaxies are $10<M_{\rm dyn}/L<1000$ and are all compatible with observed relations despite the large spread. Again, the model without feedback is an exception to the above, being very compact ($r_{1/2}=51$ pc), with a high velocity dispersion ($\sigma_{\star}=17.6~{}\rm{km}\,\rm{s}^{-1}$) and self-gravitating ($M_{\rm dyn}/L=2.0$).

We conclude that it is difficult to rule out – based solely on the observed scaling relations of dwarf galaxies – any of our models in which, at minimum, supernova feedback is included. This underlines that, while a large number of studies in the literature have been able to match such relations, they are relatively weak discriminators of the underlying physics governing the faintest dwarfs.

Given the close match to observations, it is interesting to understand whether galaxies evolve ‘along’ observed scaling relations, or if their evolutionary paths are more complex. To this end, Figure 5 shows evolutionary tracks for the quantities in the representative ‘Fiducial’ model. Dashed lines show $z=13-4$ (i.e. until star formation stops) and solid lines trace the evolution until $z=0$. The stellar velocity dispersion shows little evolution, varying by at most a factor of 3 over cosmic time. In contrast, the half mass radius evolves significantly, from $r_{1/2}<10{\,\rm pc}$ at $z>10$ to an over one order of magnitude increase in size at the current epoch. After $z=4$, most gas is completely removed due to the background UV radiation, and the galaxy expands by factor of 3, from $r_{1/2}\sim 100$ pc to $r_{1/2}\sim 300$ pc. Overall, global quantities are found to more or less evolve along the observed $z=0$ scaling relations, with the possible exception of $r_{1/2}$ and the the early evolution ($z>10$) of the estimated mass-to-light ratio.

3.4.1 The stellar mass-metallicity relation

The final scaling relation we turn to is the stellar mass (or magnitude) vs. metallicity relation (MZR). We henceforth refer to the iron to hydrogen abundance ratio ([Fe/H]) as ’metallicity’. Figure 6 shows the mean stellar [Fe/H] as a function of $M_{V}$ (left panel) and $M_{\star}$ (right panel) for all simulations, using the same point styles as in Figure 5. In the simulations, chemical abundances are calculated for each star particle following (see also Escala et al., 2018),

[TABLE]

where $Y$ and $X$ are chemical species, $m_{Y}$ and $m_{X}$ are their respective atomic masses, and $f_{X}$ and $f_{Y}$ are their respective metal mass fractions. Abundances relative to solar ($\epsilon_{X,\odot}$ and $\epsilon_{Y,\odot}$) are taken from Asplund et al. (2009) (see their table 1). The mean galactic stellar metallicity is computed, motivated by observational measurements of local group dwarf galaxies (e.g. Kirby et al., 2013), according to

[TABLE]

where [Fe/H]i is the metallicity of a star particle, $m_{\star,i}$ is the mass of a star particle and $N$ is the number of star particles in the galaxy.

The left panel contrasts our simulations with data for MW and M31 dwarf spheroidals and local group dwarf irregulars taken from Kirby et al. (2013), and UFDs from the recent review by Simon (2019)121212With UFDs defined, following the approximate magnitude separation by Simon (2019) (see his figure 5), as galaxies with $M_{V}>-7.7$.. Observations suggest a metallicity plateau around [Fe/H]$\sim-2.5$, and we indicate the lower limit observed in mean [Fe/H] (Tucana-II, [Fe/H]=$-2.90^{+0.15}_{-0.16}$, Chiti et al. 2018) with a gray dashed line in both panels. This may pose a significant challenge to current numerical simulations which tend to predict near-primordial abundances for the objects with stellar masses significantly below $10^{5}\,\mathrm{M}_{\odot}$ (right panel), or may indicate that many galaxies at these low luminositites – that are all satellites of the Milky Way and/or M31 – are tidally stripped remnants of once larger systems. We will return to this idea in future work.

By way of contrast, some of our feedback setups produce metallicities that are well-matched to observations of low mass dwarfs and UFDs, with [Fe/H]$\sim-2.7$ to $-2$ and $M_{V}\sim-6.5$ to $-5.5$. The MZR is revealed to be a highly sensitive probe of the physics related to the star formation – outflow cycle. For example, boosting the amount of energy injected by SNe leads to strong suppression of [Fe/H]. Likewise, weaker stellar feedback models lead to an overestimation of [Fe/H]. This can be understood in terms of our discussion in Section 3.3: explosive suppression of star formation using supernovae leads to expulsion of enriched gas from the interstellar medium, whereas continuous suppression of star formation, using mechanisms like UV heating, enables metals to accumulate within the interstellar medium.

It is also interesting to study the evolutionary tracks in Figure 6. No version of our galaxy evolves ‘along’ the observed relation, so that star formation quenching at any instance during these formation histories changes the match between theory and observations. Highly explosive suppression of star formation, such as that exhibited by ‘Fiducial, $10\times E_{\rm SN}$’ and ‘Fiducial, $100\times E_{\rm SN}$’, results in average metallicities131313Allowing for our Pop III metal floor to have solar abundance composition would push up these values to [Fe/H]$\sim-3$ (higher than canonical Pop III to Pop II star formation metallicity thresholds, $Z\sim 10^{-4}~{}Z_{\odot}$, Jaacks et al., 2018), which is closer to, but still in tension with, observations. We note that tests (not shown) with pop III floors at $Z\leq 10^{-4}~{}Z_{\odot}$ yield close to identical results to those presented here. This means the floor itself, at least in this simulation setup, can be introduced a posteriori to understand features such as the observed plateau in [Fe/H]. as low as [Fe/H]$\sim-4$. The opposite is true for the models with ‘weak feedback’, which rapidly evolve to almost one dex above the observed relation.

The power of the observed MZR for distinguishing the feedback models in our simulations is encouraging, since unlike the scaling relations shown in Figure 5, it may be possible to apply it as a powerful observational diagnostic of galaxy formation physics in the smallest galaxies. In future work (Orkney et al. in prep.), we will explore this using a larger simulation suite. While SN yields (here taken from Woosley & Weaver, 1995), as well as SNIa rates, will introduce uncertainties of factors $\sim 2$ at the low metallicities relevant for UFDs (Wiersma et al., 2009), these are likely sub-dominant compared to the differences between feedback implementations that we find here. A full study will also need to take into account the proximity bias of observed dwarf galaxies with sufficient spectra to derive accurate metallicities (Kirby et al., 2013) and environmental effects such as quenching due to infall onto a larger host galaxy. We discuss this further in Section 4. Finally, we note that the MZR in the faintest regime $(M_{V}\lesssim-3)$ is to date sampled by only a handful of galaxies, and in some cases mean metallicities are poorly determined (for example, Willman 1 has measured metallicities for only 2 stars). Future deep photometric observations and associated spectroscopic followups will determine the degree to which the MZR is an arbiter of galaxy formation physics.

4 Discussion

4.1 Comparison with previous work in the literature

In this section, we discuss how our results compare to recent dwarf galaxy formation results in the literature, focusing on dwarf galaxy scaling relations and the $M_{\star}-M_{200}$ relation. Only cosmological simulations reaching a spatial resolution of tens of parsecs are able to resolve the half-mass radii relevant for galaxies forming in haloes of mass $M_{200}\sim 10^{8}-10^{10}{\,\rm M}_{\odot}$ ($r_{1/2}\sim$ few 100 pc). As a result, we limit our discussion to zoom simulations (although see Rosdahl et al., 2018, for high resolution simulations of a uniform volume to $z=6$).

4.1.1 Effects of feedback on simulated dwarf galaxies

Several authors, using a variety of numerical methods, star formation and feedback prescriptions, report on global properties of simulated dwarfs forming in $M_{200}\sim 10^{9}{\,\rm M}_{\odot}$ haloes (including Wheeler et al., 2015; Oñorbe et al., 2015; Fitts et al., 2017; Macciò et al., 2017; Munshi et al., 2017; Revaz & Jablonka, 2018). The general picture is that modern simulations produce galaxies compatible with observed scaling relations, with $r_{1/2}\sim 0.1-1~{}{\,\rm kpc}$, $\sigma_{\star}\sim 5-10~{}\rm{km}\,\rm{s}^{-1}$ and stellar masses $M_{\star}\sim 10^{4}-10^{6}{\,\rm M}_{\odot}$. All successful models include subgrid models for stellar feedback that are strong enough to regulate star formation. This keeps dynamical mass-to-light ratios high, allowing stellar velocity dispersions and half-mass radii to be set by the dark matter halo alone.

While most simulations agree on this basic result, Smith et al. (2019) give an important caveat. Despite the inclusion of a momentum-capturing SN model similar to ours (see Section 2.3), their simulations all suffer from over-cooling, leading to dwarf galaxy properties similar to our ‘Fiducial, No Feedback’ case. This issue was only mitigated when stars in the Smith et al. (2019) simulations were modelled to form at a high star formation efficiency per free-fall time, $\epsilon_{\rm ff}=100\%$, leading to a more clustered injection of stellar feedback. Nonetheless, while dynamical scaling relations do capture some important aspects of the interaction between feedback algorithms (e.g. ‘blastwave’ vs. momentum based feedback) and the clustering of star formation (see e.g. the discussion in Agertz & Kravtsov 2015), we find that metal enrichment is a considerably more sensitive probe of the different feedback models that we have studied in this work.

Wheeler et al. (2015) and Oñorbe et al. (2015) used the FIRE model of feedback (Hopkins et al., 2014, see also the updated FIRE2 model, Hopkins et al. 2018), which treats a range of feedback processes, similar to our model, including a subgrid model of radiative feedback that includes both radiation pressure and ionising radiation. The Smoothed-Particle-Hydrodynamics (SPH) simulations by Macciò et al. (2017) and Revaz & Jablonka (2018) both adopted the ‘blastwave’ feedback model (Stinson et al., 2006b) that prohibits gas cooling around SNe explosions for an extended period of time, in order to allow for efficient coupling to the ISM. In addition, Macciò et al. (2017) made use of a phenomenological model for the effect of UV radiation, ‘Early Stellar Feedback’141414In this, 10% of the UV luminosity of a stellar population – typically $10^{52}$ ergs (10 times the canonical SN explosion energy) – is injected as thermal energy before any SN events take place (Stinson et al., 2013)., which is likely the reason for them finding lower stellar masses than Revaz & Jablonka (2018) ($M_{\star}\sim 10^{4}{\,\rm M}_{\odot}$ compared to $M_{\star}\sim 10^{5}{\,\rm M}_{\odot}-10^{6}{\,\rm M}_{\odot}$ for galaxies of mass $M_{200}=10^{9}{\,\rm M}_{\odot}$). The right panel of Figure 6 shows that the [Fe/H]-$M_{\star}$ relation in Macciò et al. (2017) diverge from the observed one (Kirby et al., 2014) below $M_{\star}\sim 10^{6}{\,\rm M}_{\odot}$.

Similar results were recently found by Wheeler et al. (2018) also shown in the right panel of Figure 6, at a similar numerical resolution to our ‘Hires’ simulations; for $M_{\star}\lesssim 10^{5}~{}{\,\rm M}_{\odot}$, their galaxies have mean [Fe/H]$<-3.5$, with UFDs $M_{\star}\sim 10^{3}-10^{4}~{}{\,\rm M}_{\odot}$ never enriching at all. This possibly indicates that either star formation is shut down too early in their simulated dwarfs due to overly efficient feedback, or that the IGM around star forming UFDs, at least in the local group, was significantly more pre-enriched by Pop III stars than traditionally thought, as well as predicted by galaxy simulations with Pop III enrichment (e.g. Vandenbroucke et al. 2016, but see Jaacks et al. 2018), bringing their simulated UFDs closer to the observed [Fe/H] lower limit found for faint satellite galaxies.

4.1.2 The $M_{\star}-M_{200}$ relation

At present, there is no clear theoretical consensus on expectations for $M_{\star}$ for halo masses $M_{200}$ below $\simeq 10^{10}\,{\,\rm M}_{\odot}$. While individual simulation efforts with fixed star formation and feedback prescriptions tend to produce well-defined $M_{\star}$ – $M_{200}$ relations (e.g. Wheeler et al., 2015; Fitts et al., 2017; Munshi et al., 2017; Wheeler et al., 2018), there are major differences between groups ($\sim 2$ dex, see e.g. Garrison-Kimmel et al., 2017).

The right panel of Figure 3 collects some examples from the literature. Work using the FIRE feedback and star formation models predicts steeply decreasing galaxy formation efficiencies below $M_{200}\sim 10^{10}{\,\rm M}_{\odot}$, with $M_{\star}/M_{200}<10^{-5}$ around $M_{200}\sim 10^{9}{\,\rm M}_{\odot}$ (close to a few $\times 10^{-6}$, also found in the more recent work by Wheeler et al., 2018). As such, the FIRE simulations predict the lowest galaxy formation efficiencies151515We note that a few of the Wetzel et al. (2016) ‘Latte’ dwarfs are in line with our results considered here (together with Macciò et al., 2017), possibly due to strong feedback and the inclusion of ISM heating sources such as photo-electric heating (see Section 4.2). Our ‘Hires+RT’ model predicts $M_{\star}/M_{200}\sim 5\times 10^{-5}$, one order of magnitude higher.

The data from the cosmological zoom SPH simulations in Munshi et al. (2017) show a variety of efficiencies, with results in general agreement with our entire simulation suite. By extending their work to include molecular hydrogen based star formation, akin to our approach, Munshi et al. (2018) found that stars can only form at densities high enough to allow for gas self-shielding ($n\sim 1000~{}m_{\rm H}~{}{\rm cm}^{-3}$), in agreement with our findings for the ‘Fiducial + RT, H2 based SF’ simulation. This led to a large suppression of the number of galaxies formed below $M_{200}\sim 10^{9}{\,\rm M}_{\odot}$. We found essentially no effect when adopting a H2-based prescription compared to our fiducial approach, which likely stems from us anyway restricting star formation to dense ($n>300~{}{\rm cm}^{-3}$) and cold gas ($T<100$ K) that roughly captures the environments where molecular hydrogen formation is possible. Note that we have only studied the effect of H2 physics for one particular dwarf assembly history, and a larger suite may reveal if H2 based star formation has an impact. Furthermore, the insensitivity can also stem from us already employing a high density threshold for star formation in all simulations, or that we are not modelling low enough halo masses, which we leave for a future investigation. When coupled with RT, the impact of H2-based star formation led to a $\sim 50\%$ difference in the $z=0$ stellar masses.

There is likely to be a physical spread in $M_{\star}$ at a given $M_{200}$ resulting from the diversity of mass growth histories and the different environments experienced by dwarfs (e.g. Garrison-Kimmel et al., 2017; Fitts et al., 2017; Munshi et al., 2017; Wright et al., 2018). Satellite dwarfs experience dark matter and stellar mass loss due to tidal stripping, biasing them to higher $M_{\star}/M_{200}$ ratios (Sawala et al., 2015). However, it should be emphasised that the current scatter in results from the numerical literature is almost certainly dominated by pure feedback discrepancies rather than any systematic differences in the objects being simulated.

In summary, the diversity in $M_{\star}-M_{200}$ at ‘the edge of galaxy formation’ found in the literature reflects in part specific choices in terms of numerical methods, including flavours of subgrid feedback and star formation prescriptions, as demonstrated by the large scatter found in this paper for a single choice of initial condition. In order to interpret observations robustly, both the feedback uncertainty and history variations (see also Revaz & Jablonka, 2018) will have to be taken into account; we investigate this in a companion study (Rey et al., 2019).

4.2 Simulation limitations

Having demonstrated that radiative feedback plays a crucial role in regulating the rate of star formation in ultra-faint dwarfs, we now turn to processes we have omitted. Feedback models are still far from complete, but a census of the leading-order missing ingredients and their likely effects helps to interpret the current state-of-the-art.

First, we do not model the effect of photoelectric heating (PEH), i.e. the heating by far-UV photons from young stars as they liberate electrons from dust grains (Draine, 1978). Although this effect is sub-dominant to UV heating, Forbes et al. (2016) carried out non-cosmological simulations of dwarf galaxies and found PEH to significantly reduce dwarf galaxy stellar masses, more so than supernovae feedback. On the other hand, Hu et al. (2017) suggested that the results of Forbes et al. (2016) arose due to incorrect cooling rates in the self-shielded gas. It therefore seems likely that the non-linear coupling between SNe, RT, and dispersal of dense gas, renders PEH subdominant to SN feedback. Furthermore, as the dust-to-gas ratio decreases with metallicity (super-linearly below 0.1$Z_{\odot}$, Fisher et al., 2014), PEH will most likely be a sub-dominant effect in the low metallicity environments studied in this work. Nonetheless, further work in this area is warranted.

Second, we do not include the effect of resonantly scattered Lyman-$\alpha$ photons. Self-consistently modelling the momentum transfer from the scattering of Ly-$\alpha$ photons onto the gas is computationally challenging. Kimm et al. (2018) implemented a subgrid model of this effect in RAMSES-RT, and argued that in metal poor systems, Ly-$\alpha$ photons impart momentum comparable to SNe. In their simulations of an idealised disc embedded in $M_{200}=10^{10}{\,\rm M}_{\odot}$ halo, with supernova and RT physics similar to that adopted in this work, star formation rates were reduced by a factor of two. This effect has not yet been studied in a cosmological context. Another missing feedback effect is cosmic rays (CRs). Work by Booth et al. (2013) (see also Hanasz et al., 2013; Salem & Bryan, 2014; Chan et al., 2018) demonstrated that the pressure gradients generated by cosmic rays can lead to winds with high mass loading factors, on the order of $\sim 10$, in dwarf galaxies. A full treatment CR driven winds requires magnetohydrodynamics as well as anisotropic cosmic ray diffusion coefficients. Furthermore, the impact of CRs on galactic wind properties depend strongly on the value of the diffusion coefficients, which are only empirically constrained withinn a factor of 10 for the Milky Way’s ISM, making conclusions on the impact of CRs less robust, at least currentely, than e.g. supernova feedback.

Finally, each star particle formed in our simulations is assumed to be an SSP with a Kroupa (2001) IMF. In ultra-faint galaxies, the way in which the IMF is populated likely matters, as stochasticity can affect the local level of gas heating from SNe and hence burstiness of star formation (see recent work by Applebaum et al., 2018). Although highly uncertain, if the IMF depends on the local ISM environment by e.g. becoming top-heavy at the low metallicities relevant for UFDs (as suggested by Geha et al., 2013), this can affect galaxy formation in non-trivial ways (Prgomet et al. in prep.). Furthermore, binary star physics can provide additional sources of feedback that have not been modelled here (e.g. Jeon et al., 2015).

We leave the investigation of how the above processes impact galaxy formation in low mass haloes for future work.

5 Conclusions

We have carried out cosmological zoom simulations with coupled radiation and hydrodynamics to study the formation of ultra-faint dwarf (UFD) galaxies. We studied galaxy formation in a dark matter halos, forming in relative isolation, with $z=0$ virial masses $M_{200}=10^{9}{\,\rm M}_{\odot}$, such that the cosmic UV background quenched star formation by $z\sim 4$. The simulations reached a mass and spatial resolution of $\sim 20{\,\rm M}_{\odot}$ and $\sim 3$ parsecs. Using a single realization of the Gaussian random initial conditions for our dwarf, we investigated the sensitivity of observed galaxy properties to the adopted supernova feedback model, UV background, molecular hydrogen (H2) physics, numerical resolution and multifrequency radiative transfer (RT). Our key findings are as follows:

• •

Supernova (SN) feedback lowers galaxy masses by two orders of magnitude, from $M_{\star}\sim 10^{7}{\,\rm M}_{\odot}$ when feedback is not included, to $M_{\star}\sim 10^{5}{\,\rm M}_{\odot}$. Radiative feedback, here self-consistently modeled using multi-frequency radiative transfer (RT), lowers the galaxy formation efficiency by almost an additional order of magnitude, bringing stellar masses closer to a ${\rm few}\times 10^{4}{\,\rm M}_{\odot}$, similar to local group ultra-faint dwarfs such as Eridanus-II.

• •

In models without radiative feedback, we find that dwarf galaxy formation is regulated by vigorous outflows, with much of the interstellar medium being dispersed during starburst episodes. This picture changes when RT is introduced. Radiative feedback acts to keep more gas in a warm ($\sim 10^{4}$ K) non-star forming state, leading to less gas accretion and collapse, and, as a result, less vigorous galactic scale outflows, at all times.

• •

All of our simulations with efficient feedback are in agreement with dynamical measurements of $M_{\star}$ vs. $M_{200}$ for local group dwarf satellites and UFDs, with $M_{\star}/M_{200}\sim 10^{-4}$. Molecular hydrogen based star formation and changes to the UV background modifies final stellar masses by at most a factor of two.

• •

All simulations, even those with artificially weak feedback and $M_{\star}/M_{200}>10^{-3}$, lie within the scatter of observed scaling relations for $V$-band magnitudes, half mass radii, stellar velocity dispersions and dynamical mass-to-light ratios from local group dwarf irregulars, spheroidals and UFDs, with $r_{1/2}\sim 200-400{\,\rm pc}$ and $\sigma_{\star}\sim 5-7~{}\rm{km}\,\rm{s}^{-1}$. We find that this insensitivity arises due to structural scaling relations predominantly being set by the host dark matter halo in galaxies with high mass-to-light-ratios.

• •

We find that the stellar mass-metallicity relation, based on currently available data, can differentiate galaxy formation models, as simulations fail to match observations whenever SN feedback is artificially strong or weak. This highlights that the success of galaxy formation models depends not necessarily on how many observables they can match, but rather on whether they can match a few key observables – in this case the stellar mass-metallicity relation – that are sensitive to changes in the sub-grid physics model.

Our work demonstrates that the enrichment of dwarf galaxies is a more powerful discriminant of feedback processes than any of the other observed scaling relations. The tendency for simulations of ultra-faint dwarfs in the literature to be under-enriched may suggest that their feedback models may be excessively violent, unphysically ejecting enriched gas at high redshift, or that the IGM was pre-enriched to mucher higher levels than traditionally thought. In future work (Rey et al., 2019, Orkney et al. in prep), we will study a wider range of assembly histories and halo masses in order to understand the generality of this conclusion, as well as to probe how accretion history and environment shapes both the metallicity and dynamics of the ultra-faint dwarf population.

acknowledgments

The authors thank the referee Chia-Yu Hu for his constructive comments that improved the quality of the paper. We thank Coral Wheeler, Joshua Simon, Andrew Wetzel, Andrea Macciò, Yvez Revaz and Alyson Brooks for helpful comments.

OA acknowledges support from the Swedish Research Council (grant 2014-5791) and the Knut and Alice Wallenberg Foundation. AP is supported by the Royal Society. MR acknowledges support from the Perren Fund and the IMPACT fund. JR acknowledges support from the ORAGE project from the Agence Nationale de la Recherche under grant ANR-14-CE33-0016-03.

This work was partially enabled by support from the UCL Cosmoparticle Initiative. This work was performed in part using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. This work also used the COSMA Data Centric system at Durham University, operated by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility (www.dirac.ac.uk). This equipment was funded by a BIS National E-infrastructure capital grant ST/K00042X/1, DiRAC Operations grant ST/K003267/1 and Durham University. DiRAC is part of the National E-Infrastructure. This work was also supported by a grant from the Swiss National Supercomputing Centre (CSCS) under project ID s890. Finally, a large number of simulations for this work were performed on computational resources at LUNARC, the center for scientific and technical computing at Lund University, thanks to financial support from the Royal Physiographic Society of Lund.