We present a series of high-resolution cosmological simulations1 of galaxy formation to z = 0, spanning halo masses ~108-1013M, and stellar masses ~104-1011M. Our simulations include fully explicit treatment of themultiphase interstellarmedium and stellar feedback. The stellar feedback inputs (energy, momentum, mass, and metal fluxes) are taken directly from stellar population models. These sources of feedback, with zero adjusted parameters, reproduce the observed relation between stellar and halo mass up to Mhalo ~ 1012M. We predict weak redshift evolution in the M*-Mhalo relation, consistent with current constraints to z > 6. We find that the M*-Mhalo relation is insensitive to numerical details, but is sensitive to feedback physics. Simulations with only supernova feedback fail to reproduce observed stellar masses, particularly in dwarf and high-redshift galaxies: radiative feedback (photoheating and radiation pressure) is necessary to destroy giant molecular clouds and enable efficient coupling of later supernovae to the gas. Star formation rates (SFRs) agree well with the observed Kennicutt relation at all redshifts. The galaxy-averaged Kennicutt relation is very different from the numerically imposed law for converting gas into stars, and is determined by self-regulation via stellar feedback. Feedback reduces SFRs and produces reservoirs of gas that lead to rising latetime star formation histories, significantly different from halo accretion histories. Feedback also produces large short-time-scale variability in galactic SFRs, especially in dwarfs. These properties are not captured by common 'sub-grid' wind models.

Radiative feedback (RFB) from stars plays a key role in galaxies, but remains poorly understood. We explore this using high-resolution, multifrequency radiation-hydrodynamics (RHD) simulations from the Feedback In Realistic Environments (FIRE) project. We study ultrafaint dwarf through Milky Way mass scales, including H+He photoionization; photoelectric, Lyman Werner, Compton, and dust heating; and single+multiple scattering radiation pressure (RP). We compare distinct numerical algorithms: ray-based LEBRON (exact when optically thin) and moments-based M1 (exact when optically thick). The most important RFB channels on galaxy scales are photoionization heating and single-scattering RP: in all galaxies, most ionizing/far-UV luminosity (∼1/2 of lifetime-integrated bolometric) is absorbed. In dwarfs, the most important effect is photoionization heating from the UV background suppressing accretion. In MW-mass galaxies, metagalactic backgrounds have negligible effects; but local photoionization and single-scattering RP contribute to regulating the galactic star formation efficiency and lowering central densities. Without some RFB (or other 'rapid' FB), resolved GMCs convert too-efficiently into stars, making galaxies dominated by hyperdense, bound star clusters. This makes star formation more violent and 'bursty' when SNe explode in these hyperclustered objects: thus, including RFB 'smoothes' SFHs. These conclusions are robust to RHD methods, but M1 produces somewhat stronger effects. Like in previous FIRE simulations, IR multiple-scattering is rare (negligible in dwarfs, ∼ 10 per cent of RP in massive galaxies): absorption occurs primarily in 'normal' GMCs with AV ∼ 1.

We study the implementation of mechanical feedback from supernovae (SNe) and stellar mass loss in galaxy simulations, within the Feedback In Realistic Environments (FIRE) project. We present the FIRE-2 algorithm for coupling mechanical feedback, which can be applied to any hydrodynamics method (e.g. fixed-grid, moving-mesh, and mesh-less methods), and black hole as well as stellar feedback. This algorithm ensures manifest conservation of mass, energy, and momentum, and avoids imprinting 'preferred directions' on the ejecta. We show that it is critical to incorporate both momentum and thermal energy of mechanical ejecta in a self-consistent manner, accounting for SNe cooling radii when they are not resolved. Using idealized simulations of single SN explosions, we show that the FIRE-2 algorithm, independent of resolution, reproduces converged solutions in both energy and momentum. In contrast, common 'fully thermal' (energy-dump) or 'fully kinetic' (particle-kicking) schemes in the literature depend strongly on resolution: when applied at mass resolution ≳100M⊙, they diverge by orders of magnitude from the converged solution. In galaxy-formation simulations, this divergence leads to orders-of-magnitude differences in galaxy properties, unless those models are adjusted in a resolution-dependent way. We show that all models that individually time-resolve SNe converge to the FIRE-2 solution at sufficiently high resolution (< 100M⊙). However, in both idealized single-SN simulations and cosmological galaxy-formation simulations, the FIRE-2 algorithm converges much faster than other sub-grid models without re-tuning parameters.

One important result from recent large integral field spectrograph (IFS) surveys is that the intrinsic velocity dispersion of galaxies traced by star-forming gas increases with redshift. Massive, rotation-dominated discs are already in place at z ∼ 2, but they are dynamically hotter than spiral galaxies in the local Universe. Although several plausible mechanisms for this elevated velocity dispersion (e.g. star formation feedback, elevated gas supply, or more frequent galaxy interactions) have been proposed, the fundamental driver of the velocity dispersion enhancement at high redshift remains unclear. We investigate the origin of this kinematic evolution using a suite of cosmological simulations from the FIRE (Feedback In Realistic Environments) project. Although IFS surveys generally cover a wider range of stellar masses than in these simulations, the simulated galaxies show trends between intrinsic velocity dispersion (σ intr ), SFR, and z in agreement with observations. In both observations and simulations, galaxies on the star-forming main sequence have median σ intr values that increase from z ∼ 0 to z ∼ 1–1.5, but this increasing trend is less evident at higher redshift. In the FIRE simulations, σ intr can vary significantly on time-scales of 100 Myr. These variations closely mirror the time evolution of the SFR and gas inflow rate (M gas ). By cross-correlating pairs of σ intr, M gas, and SFR, we show that increased gas inflow leads to subsequent enhanced star formation, and enhancements in σ intr tend to temporally coincide with increases in M gas and SFR.

We report the formation of bound star clusters in a sample of high-resolution cosmological zoom-in simulations of z ≥ 5 galaxies from the Feedback In Realistic Environments project. We find that bound clusters preferentially form in high-pressure clouds with gas surface densities over 104M⊙ pc-2, where the cloud-scale star formation efficiency is near unity and young stars born in these regions are gravitationally bound at birth. These high-pressure clouds are compressed by feedback-driven winds and/or collisions of smaller clouds/gas streams in highly gas-rich, turbulent environments. The newly formed clusters follow a power-law mass function of dN/dM ~ M-2. The cluster formation efficiency is similar across galaxies with stellarmasses of~107-1010M⊙ at z ≥ 5. The age spread of cluster stars is typically a few Myr and increases with cluster mass. The metallicity dispersion of cluster members is ~0.08 dex in [Z/H] and does not depend on clustermass significantly. Our findings support the scenario that present-day old globular clusters (GCs) were formed during relatively normal star formation in high-redshift galaxies. Simulations with a stricter/looser star formation model form a factor of a few more/fewer bound clusters per stellar mass formed, while the shape of the mass function is unchanged. Simulations with a lower local star formation efficiency form more stars in bound clusters. The simulated clusters are larger than observed GCs due to finite resolution. Our simulations are among the first cosmological simulations that form bound clusters self-consistently in a wide range of high-redshift galaxies.

Observations indicate that a continuous supply of gas is needed to maintain observed star formation rates in large, disky galaxies. To fuel star formation, gas must reach the inner regions of such galaxies. Despite its crucial importance for galaxy evolution, how and where gas joins galaxies is poorly constrained observationally and is rarely explored in fully cosmological simulations. To investigate gas accretion in the vicinity of galaxies, we analyze the FIRE-2 cosmological zoom-in simulations for 4 Milky Way mass galaxies (M_halo ~ 10E12 solar masses), focusing on simulations with cosmic ray physics. We find that at z~0, gas approaches the disk with angular momentum similar to the gaseous disk edge and low radial velocities, piling-up near the edge and settling into full rotational support. Accreting gas moves predominantly parallel to the disk with small but nonzero vertical velocity components, and joins the disk largely in the outskirts as opposed to "raining" down onto the disk. Once in the disk, gas trajectories are complex, being dominated by spiral arm induced oscillations and feedback. However, time and azimuthal averages show clear but slow net radial infall with transport speeds of 1-3 km/s and net mass fluxes through the disk on the order of one solar mass per year, comparable to the star formation rates of the galaxies and decreasing towards galactic center as gas is sunk into star formation. These rates are slightly higher in simulations without cosmic rays (1-7 km/s, ~4-5 solar masses per year). We find overall consistency of our results with observational constraints and discuss prospects of future observations of gas flows in and around galaxies.

We present the first measurement of the lifetimes of giant molecular clouds (GMCs) in cosmological simulations at z = 0, using the Latte suite of FIRE-2 simulations of Milky Way (MW) mass galaxies. We track GMCs with total gas mass ≳105 M⊙ at high spatial (∼1 pc), mass (7100 M⊙), and temporal (1 Myr) resolution. Our simulated GMCs are consistent with the distribution of masses for massive GMCs in the MW and nearby galaxies. We find GMC lifetimes of 5-7 Myr, or 1-2 freefall times, on average, with less than 2 per cent of clouds living longer than 20 Myr. We find decreasing GMC lifetimes with increasing virial parameter, and weakly increasing GMC lifetimes with galactocentric radius, implying that environment affects the evolutionary cycle of GMCs. However, our GMC lifetimes show no systematic dependence on GMC mass or amount of star formation. These results are broadly consistent with inferences from the literature and provide an initial investigation into ultimately understanding the physical processes that govern GMC lifetimes in a cosmological setting.

Using a state-of-the-art cosmological simulation of merging proto-galaxies at high redshift from the FIRE project, with explicit treatments of star formation and stellar feedback in the interstellar medium, we investigate the formation of star clusters and examine one of the formation hypotheses of present-day metal-poor globular clusters. We find that frequent mergers in high-redshift proto-galaxies could provide a fertile environment to produce long-lasting bound star clusters. The violent merger event disturbs the gravitational potential and pushes a large gas mass of ≳ 105-6M⊙ collectively to high density, at which point it rapidly turns into stars before stellar feedback can stop star formation. The high dynamic range of the reported simulation is critical in realizing such dense star-forming clouds with a small dynamical time-scale, tff ≲ 3 Myr, shorter than most stellar feedback time-scales. Our simulation then allows us to trace how clusters could become virialized and tightly bound to survive for up to ~420 Myr till the end of the simulation. Because the cluster's tightly bound core was formed in one short burst, and the nearby older stars originally grouped with the cluster tend to be preferentially removed, at the end of the simulation the cluster has a small age spread.

We present models of CO(1-0) emission from Milky-Way-mass galaxies at redshift zero in the FIRE-2 cosmological zoom-in simulations. We calculate the molecular abundances by post-processing the simulations with an equilibrium chemistry solver while accounting for the effects of local sources, and determine the emergent CO(1-0) emission using a line radiative transfer code. We find that the results depend strongly on the shielding length assumed, which, in our models, sets the attenuation of the incident UV radiation field. At the resolution of these simulations, commonly used choices for the shielding length, such as the Jeans length, result in CO abundances that are too high at a given H2 abundance. We find that a model with a distribution of shielding lengths, which has a median shielding length of ~3 pc in cold gas (T < 300 K) for both CO and H2, is able to reproduce both the observed CO(1-0) luminosity and inferred CO-to-H2 conversion factor at a given star formation rate compared with observations. We suggest that this short shielding length can be thought of as a subgrid model, which controls the amount of radiation that penetrates giant molecular clouds.

We analyse the different fates of the circumgalactic medium (CGM) in FIRE-2 cosmological simulations, focusing on the redshifts z = 0.25 and 2 representative of recent surveys. Our analysis includes 21 zoom-in simulations covering the halo mass range Mh(z = 0) ∼ 1010-1012 M⊙. We analyse both where the gas ends up after first leaving the CGM (its 'proximate' fate) and its location at z = 0 (its 'ultimate' fate). Of the CGM at z = 2, about half is found in the ISM or stars of the central galaxy by z = 0 in Mh(z = 2) ∼ 5 × 1011 M⊙ haloes, but most of the CGM in lower mass haloes ends up in the intergalactic medium (IGM). This is so even though most of the CGM in Mh(z = 2) ∼ 5 × 1010 M⊙ haloes first accretes on to the central galaxy before being ejected into the IGM. On the other hand, most of the CGM mass at z = 0.25 remains in the CGM by z = 0 at all halo masses analysed. Of the CGM gas that subsequently accretes on to the central galaxy in the progenitors of Mh(z = 0) ∼ 1012 M⊙ haloes, most of it is cool (T ∼ 104 K) at z = 2 but hot (∼Tvir) at z ∼ 0.25, consistent with the expected transition from cold mode to hot mode accretion. Despite the transition in accretion mode, at both z =0.25 and 2 ≥ 80 per cent of the cool gas inMh ≥ 1011 M⊙ haloes will accrete on to a galaxy. We find that the metallicity of CGM gas is typically a poor predictor of both its proximate and ultimate fates. This is because there is in general little correlation between the origin of CGM gas and its fate owing to substantial mixing while in the CGM.