We describe the second data release (DR2) of the FIRE-2 cosmological zoom-in simulations of galaxy formation, from the Feedback In Realistic Environments (FIRE) project, available at http://flathub.flatironinstitute.org/fire. DR2 includes all snapshots for most simulations, starting at z 99, with all snapshot time spacings < 25 Myr. The Core suite – comprising 14 Milky Way-mass galaxies, 5 SMC/LMC-mass galaxies, and 4 lower-mass galaxies – includes 601 snapshots to z = 0. For the Core suite, we also release resimulations with physics variations: (1) dark-matter-only versions; (2) a modified ultraviolet background with later reionization at z = 7.8; (3) magnetohydrodynamics, anisotropic conduction, and viscosity in gas; and (4) a model for cosmic-ray injection, transport, and feedback (assuming a constant diffusion coefficient). The Massive Halo suite now includes 8 massive galaxies with 278 snapshots to z = 1. The High Redshift suite includes 34 simulations: in addition to the 22 simulations run to z = 5, we now include 12 additional simulations run to z = 7 and z = 9. Most simulations include catalogs of (sub)halos and galaxies at all available snapshots, and most Core simulations to z = 0 include full halo merger trees.

The shape and orientation of dark matter (DM) halos are sensitive to the microphysics of the DM particles, yet in many mass models, the symmetry axes of the Milky Way’s DM halo are often assumed to be aligned with the symmetry axes of the stellar disk. This is well motivated for the inner DM halo, but not for the outer halo. We use zoomed-in cosmological baryonic simulations from the Latte suite of FIRE-2 Milky Way–mass galaxies to explore the evolution of the DM halo’s orientation with radius and time, with or without a major merger with a Large Magellanic Cloud analog, and when varying the DM model. In three of the four cold DM halos we examine, the orientation of the halo minor axis diverges from the stellar disk vector by more than 20° beyond about 30 galactocentric kpc, reaching a maximum of 30°–90°, depending on the individual halo’s formation history. In identical simulations using a model of self-interacting DM with σ = 1 cm 2 g −1 , the halo remains aligned with the stellar disk out to ∼200–400 kpc. Interactions with massive satellites ( M ≳ 4 × 10 10 M ⊙ at pericenter; M ≳ 3.3 × 10 10 M ⊙ at infall) affect the orientation of the halo significantly, aligning the halo’s major axis with the satellite galaxy from the disk to the virial radius. The relative orientation of the halo and disk beyond 30 kpc is a potential diagnostic of self-interacting DM, if the effects of massive satellites can be accounted for.

We perform cosmological hydrodynamical simulations to study the formation of proto-globular cluster candidates in progenitors of present-day dwarf galaxies ( M-vir approximate to 10(10) M-circle dot at z = 0) as part of the 'Feedback in Realistic Environment' (FIRE) project. Compact ( r (1/2) < 30 pc), relatively massive (0.5 x 10 (5) less than or similar to M-*/M-circle dot less than or similar to 5 x 10 (5) ), self-bound stellar clusters form at 11 greater than or similar to 5 in progenitors with M (vi)r approximate to 10(9) M-circle dot. Cluster formation is triggered when at least 10 7 Me of dense, turbulent gas reaches Egas 104 Me pc-2 as a result of the compressive effects of supernova feedback or from cloud-cloud collisions. The clusters can survive for 2 - 3 Gyr ; absent numerical effects, they could possibly survive substantially longer, perhaps to z = 0. The longest lived clusters are those that form at significant distance - several hundreds of pc - from their host galaxy. We therefore predict that globular clusters forming in progenitors of present-day dwarf galaxies will be offset from any pre-existing stars within their host dark matter haloes as opposed to deeply embedded within a well-defined galaxy. Properties of the nascent clusters are consistent with observations of some of the faintest and most compact high-redshift sources in Hubble Space Telescope lensing fields and are at the edge of what will be detectable as point sources in deep imaging of non-lensed fields with JWST. By contrast, the star clusters' host galaxies will remain undetectable.

We describe a public data release of the FIRE-2 cosmological zoom-in simulations of galaxy formation (available at http://flathub.flatironinstitute.org/fire ) from the Feedback In Realistic Environments (FIRE) project. FIRE-2 simulations achieve parsec-scale resolution to explicitly model the multiphase interstellar medium while implementing direct models for stellar evolution and feedback, including stellar winds, core-collapse and Type Ia supernovae, radiation pressure, photoionization, and photoelectric heating. We release complete snapshots from three suites of simulations. The first comprises 20 simulations that zoom in on 14 Milky Way (MW)–mass galaxies, five SMC/LMC-mass galaxies, and four lower-mass galaxies including one ultrafaint; we release 39 snapshots across z = 0–10. The second comprises four massive galaxies, with 19 snapshots across z = 1–10. Finally, a high-redshift suite comprises 22 simulations, with 11 snapshots across z = 5–10. Each simulation also includes dozens of resolved lower-mass (satellite) galaxies in its zoom-in region. Snapshots include all stored properties for all dark matter, gas, and star particles, including 11 elemental abundances for stars and gas, and formation times (ages) of star particles. We also release accompanying (sub)halo catalogs, which include galaxy properties and member star particles. For the simulations to z = 0, including all MW-mass galaxies, we release the formation coordinates and an “ex situ” flag for all star particles, pointers to track particles across snapshots, catalogs of stellar streams, and multipole basis expansions for the halo mass distributions. We describe publicly available python packages for reading and analyzing these simulations.

Increasingly, uncertainties in predictions from galaxy formation simulations (at sub-Milky Way masses) are dominated by uncertainties in stellar evolution inputs. In this paper, we present the full set of updates from the FIRE-2 version of the Feedback In Realistic Environments (FIRE) project code, to the next version, FIRE-3. While the transition from FIRE-1 to FIRE-2 focused on improving numerical methods, here we update the stellar evolution tracks used to determine stellar feedback inputs, e.g. stellar mass-loss (O/B and AGB), spectra (luminosities and ionization rates), and supernova rates (core-collapse and Ia), as well as detailed mass-dependent yields. We also update the low-temperature cooling and chemistry, to enable improved accuracy at $T \lesssim 10^{4}\,$K and densities $n\gg 1\,{\rm cm^{-3}}$, and the meta-galactic ionizing background. All of these synthesize newer empirical constraints on these quantities and updated stellar evolution and yield models from a number of groups, addressing different aspects of stellar evolution. To make the updated models as accessible as possible, we provide fitting functions for all of the relevant updated tracks, yields, etc, in a form specifically designed so they can be directly 'plugged in' to existing galaxy formation simulations. We also summarize the default FIRE-3 implementations of 'optional' physics, including spectrally-resolved cosmic rays and supermassive black hole growth and feedback.

We study how supersonic streaming velocities of baryons relative to dark matter—a large-scale effect imprinted at recombination and coherent over ∼3 Mpc scales—affect the formation of dwarf galaxies at z ≳ 5. We perform cosmological hydrodynamic simulations, including and excluding streaming velocities, in regions centered on halos with M _vir ( z = 0) ≈ 10 ^10 M _⊙ ; the simulations are part of the Feedback In Realistic Environments (FIRE) project and run with FIRE-3 physics. Our simulations comprise many thousands of systems with halo masses between M _vir = 2 × 10 ^5 M _⊙ and 2 × 10 ^9 M _⊙ in the redshift range z = 20–5. A few hundred of these galaxies form stars and have stellar masses ranging from 100 to 10 ^7 M _⊙ . While star formation is globally delayed by approximately 50 Myr in the streaming relative to nonstreaming simulations and the number of luminous galaxies is correspondingly suppressed at high redshift in the streaming runs, these effects decay with time. By z = 5, the properties of the simulated galaxies are nearly identical in the streaming versus nonstreaming runs, indicating that any effects of streaming velocities on the properties of galaxies at the mass scale of classical dwarfs and larger do not persist to z = 0.

Self-interacting dark matter (SIDM) models offer one way to reconcile inconsistencies between observations and predictions from collisionless cold dark matter (CDM) models on dwarf-galaxy scales. In order to incorporate the effects of both baryonic and SIDM interactions, we study a suite of cosmological-baryonic simulations of Milky-Way (MW)-mass galaxies from the Feedback in Realistic Environments (FIRE-2) project where we vary the SIDM self-interaction cross-section σ/m. We compare the shape of the main dark matter (DM) halo at redshift z=0 predicted by SIDM simulations (at σ/m=0.1, 1, and 10 cm^2 g^-1) with CDM simulations using the same initial conditions. In the presence of baryonic feedback effects, we find that SIDM models do not produce the large differences in the inner structure of MW-mass galaxies predicted by SIDM-only models. However, we do find that the radius where the shape of the total mass distribution begins to differ from that of the stellar mass distribution is dependent on σ/m. This transition could potentially be used to set limits on the SIDM cross-section in the MW.

ABSTRACT Self-interacting dark matter (SIDM) models have received great attention over the past decade as solutions to the small-scale puzzles of astrophysics. Though there are different implementations of dark matter (DM) self-interactions in N-body codes of structure formation, there has not been a systematic study to compare the predictions of these different implementations. We investigate the implementation of dark matter self-interactions in two simulation codes:gizmo and arepo. We begin with identical initial conditions for an isolated 1010 M⊙ dark matter halo and investigate the evolution of the density and velocity dispersion profiles in gizmo and arepo for SIDM cross-section over mass of 1, 5, and 50 $\rm cm^2\, g^{-1}$. Our tests are restricted to the core expansion phase, where the core density decreases and core radius increases with time. We find better than 30 per cent agreement between the codes for the density profile in this phase of evolution, with the agreement improving at higher resolution. We find that varying code-specific SIDM parameters changes the central halo density by less than 10 per cent outside of the convergence radius. We argue that SIDM core formation is robust across the two different schemes and conclude that these codes can reliably differentiate between cross-sections of 1, 5, and 50 $\rm cm^2\, g^{-1}$, but finer distinctions would require further investigation.

We present a suite of baryonic cosmological zoom-in simulations of self-interacting dark matter (SIDM) haloes within the ‘Feedback In Realistic Environment’ (FIRE) project. The three simulated haloes have virial masses of $\sim 10^{12}\, \text{M}_\odot$ at z = 0, and we study velocity-independent self-interaction cross sections of 1 and 10 ${\rm cm^2 \, g^{-1}}$. We study star formation rates and the shape of dark matter density profiles of the parent haloes in both cold dark matter (CDM) and SIDM models. Galaxies formed in the SIDM haloes have higher star formation rates at z ≤ 1, resulting in more massive galaxies compared to the CDM simulations. While both CDM and SIDM simulations show diverse shape of the dark matter density profiles, the SIDM haloes can reach higher and more steep central densities within few kpcs compared to the CDM haloes. We identify a correlation between the build-up of the stars within the half-mass radii of the galaxies and the growth in the central dark matter densities. The thermalization process in the SIDM haloes is enhanced in the presence of a dense stellar component. Hence, SIDM haloes with highly concentrated baryonic profiles are predicted to have higher central dark matter densities than the CDM haloes. Overall, the SIDM haloes are more responsive to the presence of a massive baryonic distribution than their CDM counterparts.

Core formation and runaway core collapse in models with self-interacting dark matter (SIDM) significantly alter the central density profiles of collapsed halos. Using a forward modeling inference framework with simulated datasets, we demonstrate that flux ratios in quadruple image strong gravitational lenses can detect the unique structural properties of SIDM halos, and statistically constrain the amplitude and velocity dependence of the interaction cross section in halos with masses between $10^6 - 10^{10} M_{\odot}$. Measurements on these scales probe self-interactions at velocities below $30 \ \rm{km} \ \rm{s^{-1}}$, a relatively unexplored regime of parameter space, complimenting constraints at higher velocities from galaxies and clusters. We cast constraints on the amplitude and velocity dependence of the interaction cross section in terms of $\sigma_{20}$, the cross section amplitude at $20 \ \rm{km} \ \rm{s^{-1}}$. With 50 lenses, a sample size available in the near future, and flux ratios measured from spatially compact mid-IR emission around the background quasar, we forecast $\sigma_{20} < 11-23 \ \rm{cm^2} \rm{g^{-1}}$ at $95 \%$ CI, depending on the amplitude of the subhalo mass function, and assuming cold dark matter (CDM). Alternatively, if $\sigma_{20} = 19.2 \ \rm{cm^2}\rm{g^{-1}}$ we can rule out CDM with a likelihood ratio of 20:1, assuming an amplitude of the subhalo mass function that results from doubly-efficient tidal disruption in the Milky Way relative to massive elliptical galaxies. These results demonstrate that strong lensing of compact, unresolved sources can constrain SIDM structure on sub-galactic scales across cosmological distances, and the evolution of SIDM density profiles over several Gyr of cosmic time.

We study evolution of self-interacting dark matter subhalos in the Milky Way tidal field. The interaction between the subhalos and the Milky Way's tides lead to more diverse dark matter distributions in the inner region, compared to their cold dark matter counterparts. We test this scenario with two Milky Way satellite galaxies, Draco and Fornax, opposite extremes in the inner dark matter content, and find that they can be accommodated within the self-interacting dark matter model proposed to explain the diverse rotation curves of spiral galaxies in the field.

We perform a series of controlled N-body simulations to study realizations of the recently discovered Antlia 2 galaxy in cold dark matter (CDM) and self-interacting dark matter (SIDM) scenarios. Our simulations contain six benchmark models, where we vary the initial halo concentration and the self-scattering cross section. We adopt well-motivated initial stellar and halo masses, and our fiducial orbit has a small pericenter. After evolving in the Milky Way's tidal field, the simulated galaxies experience significant mass loss and their stellar distributions expand accordingly. These tidal effects are more prominent if the initial halo concentration is lower and if the self-scattering cross section is larger. Our results show that Antlia 2-like galaxies could be realized in CDM if the halo concentration is low and the stellar distribution is diffuse at the infall time, while these conditions could be relaxed in SIDM. We also find all the simulated galaxies predict approximately the same stellar velocity dispersion after imposing selection criteria for stellar particles. This has important implications for testing dark matter models using tidally disturbed systems.

We explore a scenario where metal poor globular clusters (GCs) form at the centres of their own dark matter halos in the early universe before reionization. This hypothesis leads to predictions about the abundance, distribution and kinematics of GCs today that we explore using cosmological N-body simulations and analytical modelling. We find that selecting the massive tail of collapsed objects at $z\gtrsim 9$ as GC formation sites leads to four main predictions: i) a highly clustered population of GCs around galaxies today, ii) a natural scaling between number of GCs and halo virial mass that follows roughly the observed trend, iii) a very low number of free floating GCs outside massive halos and iv) GCs should be embedded within massive and extended dark matter (sub)halos. We find that the strongest constraint to the model is given by the combination of (i) and (ii): a mass cut to tagged GCs halos which accounts for the number density of metal poor GCs today predicts a radial distribution that is too extended compared to recent observations. On the other hand, a mass cut sufficient to match the observed half number radius could only explain 60% of the metal poor population. In all cases, observations favour early redshifts for GC formation ($z\geq 15$).

We study halo mass functions with the Press–Schechter formalism for interacting dark matter models, where matter power spectra are damped due to dark acoustic oscillations in the early universe. After adopting a smooth window function, we calibrate the analytical model with numerical simulations from the “effective theory of structure formation” (ETHOS) project and fix the model parameters in the high-mass regime, Mh ≳ 3 × 1010 M⊙. We also perform high-resolution cosmological simulations with halo masses down to Mh ∼ 108 M⊙ to cover a wide mass range for comparison. Although the model is calibrated with ETHOS1 and cold dark matter simulations for high halo masses at redshift z = 0, it successfully reproduces simulations for two other ETHOS models in the low-mass regime at low and high redshifts. As an application, we compare the cumulative number density of halos to that of observed galaxies at z = 6, and find that the interacting dark matter models with a kinetic decoupling temperature below 0.5 keV is disfavored. We also perform the abundance-matching analysis and derive the stellar–halo mass relation for these models at z = 4. Suppression in halo abundance leads to less massive halos that host observed galaxies in the stellar mass range M * ≃ 10 5 – 10 7 M ⊙ .

We employ isolated N-body simulations to study the response of self-interacting dark matter (SIDM) halos in the presence of the baryonic potentials. Dark matter self-interactions lead to kinematic thermalization in the inner halo, resulting in a tight correlation between the dark matter and baryon distributions. A deep baryonic potential shortens the phase of SIDM core expansion and triggers core contraction. This effect can be further enhanced by a large self-scattering cross section. We find the final SIDM density profile is sensitive to the baryonic concentration and the strength of dark matter self-interactions. Assuming a spherical initial halo, we also study evolution of the SIDM halo shape together with the density profile. The halo shape at later epochs deviates from spherical symmetry due to the influence of the non-spherical disc potential, and its significance depends on the baryonic contribution to the total gravitational potential, relative to the dark matter one. In addition, we construct a multi-component model for the Milky Way, including an SIDM halo, a stellar disc and a bulge, and show it is consistent with observations from stellar kinematics and streams.

Galactic rotation curves are a fundamental constraint for any cosmological model. We use controlled N-body simulations of galaxies to study the gravitational effect of baryons in a scenario with collisionless cold dark matter (CDM) versus one with a self-interacting dark matter (SIDM) component. In particular, we examine the inner profiles of the rotation curves in the velocity range Vmax = 30-250 km/s, whose diversity has been found to be greater than predicted by the Lambda-CDM scenario. We find that the scatter in the observed rotation curves exceeds that predicted by dark matter only mass-concentration relations in either the CDM nor SIDM models. Allowing for realistic baryonic content and spatial distributions, however, helps create a large variety of rotation curve shapes; which is in better agreement with observations in the case of self-interactions due to the characteristic cored profiles being more accommodating to the slowly rising rotation curves than CDM. We find individual fits to model two of the most remarkable outliers of similar Vmax, UGC 5721 and IC 2574, the former a cusp-like rotation curve and the latter a seemingly 8 kpc cored profile. This diversity in SIDM arises as permutations of overly concentrated halos with compact baryonic distributions versus underdense halos with extended baryonic disks. The SIDM solution is promising and its feasibility ultimately depends on the sampling of the halo mass-concentration relation and its interplay with the baryonic profiles, emphasising the need for a better understanding of the frequency of extreme outliers present in current observational samples.