Cosmological Zoom-in Simulations Research Articles

Black hole–galaxy scaling relations in FIRE: the importance of black hole location and mergers

ABSTRACT The concurrent growth of supermassive black holes (SMBHs) and their host galaxies remains to be fully explored, especially at high redshift. While often understood as a consequence of self-regulation via AGN feedback, it can also be explained by alternative SMBH accretion models. Here, we expand on previous work by studying the growth of SMBHs with the help of a large suite of cosmological zoom-in simulations (MassiveFIRE) that are part of the Feedback in Realistic Environments (FIRE) project. The growth of SMBHs is modelled in post-processing with different black hole accretion models, placements, and merger treatments, and validated by comparing to on-the-fly calculations. Scaling relations predicted by the gravitational torque-driven accretion (GTDA) model agree with observations at low redshift without the need for AGN feedback, in contrast to models in which the accretion rate depends strongly on SMBH mass. At high redshift, we find deviations from the local scaling relations in line with previous theoretical results. In particular, SMBHs are undermassive, presumably due to stellar feedback, but start to grow efficiently once their host galaxies reach M* ∼ 1010M⊙. We analyse and explain these findings in the context of a simple analytic model. Finally, we show that the predicted scaling relations depend sensitively on the SMBH location and the efficiency of SMBH merging, particularly in low-mass systems. These findings highlight the relevance of understanding the evolution of SMBH-galaxy scaling relations to predict the rate of gravitational wave signals from SMBH mergers across cosmic history.

High-redshift quasars and their host galaxies – II. Multiphase gas and stellar kinematics

ABSTRACT Observations of z ≳ 6 quasars provide information on the early phases of the most massive black holes (MBHs) and galaxies. Current observations at sub-mm wavelengths trace cold and warm gases, and future observations will extend information to other gas phases and the stellar properties. The goal of this study is to examine the gas life cycle in a z ≳ 6 quasar: from accretion from the halo to the galaxy and all the way into the MBH, to how star formation and the MBH itself affect the gas properties. Using a very high resolution cosmological zoom-in simulation of a z = 7 quasar, including state-of-the-art non-equilibrium chemistry, MBH formation, growth, and feedback, we investigate the distribution of the different gas phases in the interstellar medium across cosmic time. We assess the morphological evolution of the quasar host using different tracers (star- or gas-based) and the thermodynamic distribution of the MBH accretion-driven outflows, finding that obscuration in the disc is mainly due to molecular gas, with the atomic component contributing at larger scales and/or above/below the disc plane. Moreover, our results also show that molecular outflows, if present, are more likely the result of gas being lifted near the MBH than production within the wind because of thermal instabilities. Finally, we also discuss how different gas phases can be employed to dynamically constrain the MBH mass, and argue that resolutions below ∼100 pc yield unreliable estimates because of the strong contribution of the nuclear stellar component to the potential at larger scales.

Gas Dynamics in the Galaxy: Total Mass Distribution and the Bar Pattern Speed

Abstract Gas morphology and kinematics in the Milky Way contain key information for understanding the formation and evolution of our Galaxy. We present hydrodynamical simulations based on realistic barred Milky Way potentials constrained by recent observations. Our model can reproduce most features in the observed longitude–velocity diagram, including the Central Molecular Zone, the Near and Far 3 kpc arms, the Molecular Ring, and the spiral arm tangents. It can also explain the noncircular motions of masers from the recent BeSSeL2 survey. The central gas kinematics are consistent with a mass of 6.9 × 108 M ⊙ in the Nuclear Stellar Disk. Our model predicts the formation of an elliptical gaseous ring surrounding the bar, which is composed of the 3 kpc arms, the Norma arm, and the bar-spiral interfaces. This ring is similar to those “inner” rings in some Milky Way analogs with a boxy/peanut-shaped bulge (e.g., NGC 4565 and NGC 5746). The kinematics of gas near the solar neighborhood are governed by the Local arm. The bar pattern speed constrained by our gas model is 37.5–40 km s−1 kpc−1, corresponding to a corotation radius of R CR = 6.0–6.4 kpc. The rotation curve of our model rises gently within the central ∼ 5 kpc, significantly less steep than those predicted by some recent zoom-in cosmological simulations.

Simulating Globular Clusters in Dark Matter Subhalos

Abstract A cosmological zoom-in simulation that develops into a Milky Way-like halo begins at redshift 7. The initial dark matter distribution is seeded with dense star clusters of median mass 5 × 105 M ⊙, placed in the largest subhalos present, which have a median peak circular velocity of 25 km s−1. Three simulations are initialized using the same dark matter distribution with the star clusters starting on approximately circular orbits having initial median radii 6.8, 0.14 kpc, and, at the exact center of the subhalos. The simulations are evolved to the current epoch at which time the median galactic orbital radii of the three sets of clusters are 30, 5, and 16 kpc, with the clusters losing about 2%, 50%, and 15% of their mass, respectively. Clusters beginning at small orbital radii have so much tidal forcing that they are often not in equilibrium. Clusters that start at larger subhalo radii have a velocity dispersion that declines smoothly to ≃20% of the central value at ≃20 half-mass radii. The clusters that begin in the subhalo centers can show a rise in velocity dispersion beyond 3–5 half-mass radii. That is, the clusters that form without local dark matter always have stellar-mass-dominated kinematics at all radii, whereas about 25% of the clusters that begin in subhalo centers have remnant local dark matter.

On the origin of surprisingly cold gas discs in galaxies at high redshift

ABSTRACT We address the puzzling observational indications for very ‘cold’ galactic discs at redshifts z ≳ 3, an epoch when discs are expected to be highly perturbed. Using a high-resolution cosmological zoom-in simulation, we identify such a cold disc at z ∼ 3.5, with a rotation velocity to velocity dispersion ratio of vϕ/σr ≃ 5 for the total gas. It forms as a result of a period of intense accretion of co-planar, co-rotating gas via cold cosmic-web streams. This thin disc survives for ∼5 orbital periods, after which it is disrupted by mergers and counter-rotating streams, longer but consistent with our estimate that a galaxy of this mass ($M_\star \sim 10^{10}\, \mathrm{M_\odot }$) typically survives merger-driven spin flips for ∼2–3 orbital periods. We find that vϕ/σr is highly sensitive to the tracer used to perform the kinematic analysis. While it is vϕ/σr ≃ 3.5 for atomic H i gas, it is vϕ/σr ≃ 8 for molecular CO and H2. This reflects the confinement of molecular gas to cold, dense clouds that reside near the disc mid-plane, while the atomic gas is spread into a turbulent and more extended thicker disc. The proposed mechanisms is a theoretical proposal that has not been validated yet with proper statistical measurements and it remains unclear whether it occurs frequently enough to explain the multiple discoveries of cold gas discs in high-z galaxies.

Impact of gas spin and Lyman–Werner flux on black hole seed formation in cosmological simulations: implications for direct collapse

ABSTRACT Direct collapse black holes (BHs) are promising candidates for producing massive z ≳ 6 quasars, but their formation requires fine-tuned conditions. In this work, we use cosmological zoom simulations to study systematically the impact of requiring: (1) low gas angular momentum (spin), and (2) a minimum incident Lyman–Werner (LW) flux in order to form BH seeds. We probe the formation of seeds (with initial masses of $M_{\rm seed} \sim 10^4\!-\!10^6\, \mathrm{M}_{\odot }\, h^{-1})$ in haloes with a total mass >3000 × Mseed and a dense, metal-poor gas mass >5 × Mseed. Within this framework, we find that the seed-forming haloes have a prior history of star formation and metal enrichment, but they also contain pockets of dense, metal-poor gas. When seeding is further restricted to haloes with low gas spins, the number of seeds formed is suppressed by factors of ∼6 compared to the baseline model, regardless of the seed mass. Seed formation is much more strongly impacted if the dense, metal-poor gas is required to have a critical LW flux (Jcrit). Even for Jcrit values as low as 50J21, no $8\times 10^{5}~\mathrm{M}_{\odot }\, h^{-1}$ seeds are formed. While lower mass ($1.25\times 10^{4},1\times 10^{5}~\mathrm{M}_{\odot }\, h^{-1}$) seeds do form, they are strongly suppressed (by factors of ∼10–100) compared to the baseline model at gas mass resolutions of $\sim 10^4~\mathrm{M}_{\odot }\, h^{-1}$ (with even stronger suppression at higher resolutions). As a result, BH merger rates are also similarly suppressed. Since early BH growth is dominated by mergers in our models, none of the seeds are able to grow to the supermassive regime ($\gtrsim 10^6~\mathrm{M}_{\odot }\, h^{-1}$) by z = 7. Our results hint that producing the bulk of the z ≳ 6 supermassive BH population may require alternate seeding scenarios that do not depend on the LW flux, early BH growth dominated by rapid or super-Eddington accretion, or a combination of these possibilities.

Testing for Dark Matter in the Outskirts of Globular Clusters

Abstract The proper motions of stars in the outskirts of globular clusters are used to estimate cluster velocity dispersion profiles as far as possible within their tidal radii. We use individual color–magnitude diagrams to select high-probability cluster stars for 25 metal-poor globular clusters within 20 kpc of the Sun, 19 of which have substantial numbers of stars at large radii. Of the 19, 11 clusters have a falling velocity dispersion in the 3–6 half-mass radii range, 6 are flat, and 2 plausibly have a rising velocity dispersion. The profiles are all in the range expected from simulated clusters that started at high redshift in a zoom-in cosmological simulation. The 11 clusters with falling velocity dispersion profiles are consistent with no dark matter above the Galactic background. The six clusters with approximately flat velocity dispersion profiles could have local dark matter, but are ambiguous. The two clusters with rising velocity dispersion profiles are consistent with a remnant local dark matter halo, but need membership confirmation and detailed orbital modeling to further test these preliminary results.

Too dense to go through: the role of low-mass clusters in the pre-processing of satellite galaxies

ABSTRACT We study the evolution of satellite galaxies in clusters of the c-eagle simulations, a suite of 30 high-resolution cosmological hydrodynamical zoom-in simulations based on the eagle code. We find that the majority of galaxies that are quenched at z = 0 (≳80${{\ \rm per\ cent}}$) reached this state in a dense environment (log10M200[M⊙] ≥13.5). At low redshift, regardless of the final cluster mass, galaxies appear to reach their quenching state in low-mass clusters. Moreover, galaxies quenched inside the cluster that they reside in at z = 0 are the dominant population in low-mass clusters, while galaxies quenched in a different halo dominate in the most massive clusters. When looking at clusters at z > 0.5, their in situ quenched population dominates at all cluster masses. This suggests that galaxies are quenched inside the first cluster they fall into. After galaxies cross the cluster’s r200 they rapidly become quenched (≲1 Gyr). Just a small fraction of galaxies ($\lesssim 15{{\ \rm per\ cent}}$) is capable of retaining their gas for a longer period of time, but after 4 Gyr, almost all galaxies are quenched. This phenomenon is related to ram pressure stripping and is produced when the density of the intracluster medium reaches a threshold of $\rho _{\rm ICM}\, \sim 3 \times 10 ^{-5}$ nH (cm−3). These results suggest that galaxies start a rapid-quenching phase shortly after their first infall inside r200 and that, by the time they reach r500, most of them are already quenched.

Gas infall and radial transport in cosmological simulations of milky way-mass discs

ABSTRACT Observations indicate that a continuous supply of gas is needed to maintain observed star formation rates in large, discy 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 rarely explored in fully cosmological simulations. To investigate gas accretion in the vicinity of galaxies at low redshift, we analyse the FIRE-2 cosmological zoom-in simulations for 4 Milky Way mass galaxies (Mhalo ∼ 1012M⊙), focusing on simulations with cosmic ray physics. We find that at z ∼ 0, gas approaches the disc with angular momentum similar to the gaseous disc edge and low radial velocities, piling-up near the edge and settling into full rotational support. Accreting gas moves predominately parallel to the disc and joins largely in the outskirts. Immediately prior to joining the disc, trajectories briefly become more vertical on average. Within the disc, gas motion is complex, being dominated by spiral arm induced oscillations and feedback. However, time and azimuthal averages show slow net radial infall with transport speeds of 1–3 km s−1 and net mass fluxes through the disc of ∼M⊙ yr−1, comparable to the galaxies’ star formation rates and decreasing towards galactic centre as gas is sunk into star formation. These rates are slightly higher in simulations without cosmic rays (1–7 km s−1, ∼4–5 M⊙ yr−1). We find overall consistency of our results with observational constraints and discuss prospects of future observations of gas flows in and around galaxies.

Gas flows in galaxy mergers: supersonic turbulence in bridges, accretion from the circumgalactic medium, and metallicity dilution

ABSTRACT In major galaxy mergers, the orbits of stars are violently perturbed, and gas is torqued to the centre, diluting the gas metallicity and igniting a starburst. In this paper, we study the gas dynamics in and around merging galaxies using a series of cosmological magnetohydrodynamical zoom-in simulations. We find that the gas bridge connecting the merging galaxies pre-coalescence is dominated by turbulent pressure, with turbulent Mach numbers peaking at values of 1.6–3.3. This implies that bridges are dominated by supersonic turbulence, and are thus ideal candidates for studying the impact of extreme environments on star formation. We also find that gas accreted from the circumgalactic medium (CGM) during the merger significantly contributes (27–51 per cent) to the star formation rate (SFR) at the time of coalescence and drives the subsequent reignition of star formation in the merger remnant. Indeed, 19–53 per cent of the SFR at $z$ = 0 originates from gas belonging to the CGM prior the merger. Finally, we investigate the origin of the metallicity-diluted gas at the centre of merging galaxies. We show that this gas is rapidly accreted on to the Galactic Centre with a time-scale much shorter than that of normal star-forming galaxies. This explains why coalescing galaxies are not well-captured by the fundamental metallicity relation.

Panspermia in a Milky Way–like Galaxy

Abstract We study the process of panspermia in Milky Way–like galaxies by modeling the probability of successful travel of organic compounds between stars harboring potentially habitable planets. To this end, we apply the modified habitability recipe of Gobat & Hong to a model galaxy from the McMaster Unbiased Galaxy Simulations suite of zoom-in cosmological simulations. We find that, unlike habitability, which only occupies a narrow dynamic range over the entire galaxy, the panspermia probability can vary by orders of magnitude between the inner (R, b = 1–4 kpc) and outer disk. However, only a small fraction of star particles have very large values for the panspermia probability and, consequently, the fraction of star particles where the panspermia process is more effective than prebiotic evolution is much lower than from naïve expectations based on the ratio between the panspermia probability and natural habitability.

Population III star formation in an X-ray background – I. Critical halo mass of formation and total mass in stars

ABSTRACT The first luminous objects forming in the universe produce radiation backgrounds in the far-ultraviolet and X-ray bands that affect the formation of Population III stars. Using a grid of cosmological hydrodynamics zoom-in simulations, we explore the impact of the Lyman–Warner (LW) and X-ray radiation backgrounds on the critical dark matter (DM) halo mass for Population III star formation and the total mass in stars per halo. We find that the LW radiation background lowers the H2 fraction and delays the formation of the Population III stars. On the other hand, X-ray irradiation anticipates the redshift of collapse and reduces the critical halo mass, unless the X-ray background is too strong and gas heating shuts down gas collapse into the haloes and prevents star formation. Therefore, an X-ray background can increase the number of DM haloes forming Population III stars by about a factor of 10, but the total mass in stars forming in each halo is reduced. This is because X-ray radiation increases the molecular fraction and lowers the minimum temperature of the collapsing gas (or equivalently the mass of the quasi-hydrostatic core) and therefore slows down the accretion of the gas on to the central protostar.

The Effects of Dark Matter and Baryonic Physics on the Milky Way Subhalo Population in the Presence of the Large Magellanic Cloud

Abstract Given recent developments in our understanding of the Large Magellanic Cloud's (LMC) impact on the Milky Way's (MW) dark matter subhalo population, we compare the signatures of dark matter and baryonic physics on subhalos in MW systems with realistic LMC analogs. In particular, we study the effects of self-interacting dark matter (SIDM), warm dark matter (WDM), and the Galactic disk on the peak maximum circular velocity (V peak) function, radial distribution, and spatial distribution of MW and LMC-associated subhalos using cosmological dark matter-only zoom-in simulations of MW+LMC systems. For a fixed abundance of subhalos expected to host dwarf galaxies (V peak ≳ 20 km s−1), SIDM and WDM can produce a similar mass-dependent suppression of the subhalo V peak function, while disk disruption is mass independent. Subhalos in the inner regions of the MW are preferentially disrupted by both self-interactions and the disk, while suppression in WDM is radially independent. The relative abundance of LMC-associated subhalos is not strongly affected by disk disruption or WDM, but is significantly suppressed in SIDM due to self-interactions with the LMC at early times and with the MW during LMC infall at late times, erasing spatial anisotropy in the MW subhalo population. These results provide avenues to distinguish dark matter and baryonic physics by combining properties of the MW and LMC subhalo populations probed by upcoming observations of satellite galaxies and stellar streams.

Turbulence in the intragroup and circumgalactic medium

Context. In massive objects, such as galaxy clusters, the turbulent velocity dispersion, σturb, is tightly correlated to both the object mass, M, and the thermal energy. Aims. Here, we investigate whether these scaling laws extend to lower-mass objects in dark-matter filaments. Methods. We perform a cosmological zoom-in simulation of a filament using an adaptive filtering technique for the resolved velocity component and a subgrid-scale model to account for the unresolved component. We then compute the mean turbulent and thermal energies for all halos in the zoom-in region and compare different definitions of halo averages. Averaging constrained by density and temperature thresholds is favored over averages solely based on virial spheres. Results. We find no clear trend for the turbulent velocity dispersion versus halo mass, but significant correlation and a scaling law with exponent α ∼ 0.5 between the turbulent velocity dispersion and thermal energy that agrees with a nearly constant turbulent Mach number, similar to more massive objects. Conclusions. We conclude that the self-similar energetics proposed for galaxy clusters extends down to the circumgalactic medium of individual galaxies.

Some First Stars Were Red: Detecting Signatures of Massive Population III Formation through Long-term Stochastic Color Variations

Abstract Identifying stars formed in pristine environments (Population III) within the first billion years is vital to uncovering the earliest growth and chemical evolution of galaxies. Population III galaxies, however, are typically expected to be too faint and too few in number to be detectable by forthcoming instruments without extremely long integration times and/or extreme lensing. In an environment, however, where star formation is suppressed until a halo crosses the atomic-cooling limit (e.g., by a modest Lyman–Werner flux, high baryonic streaming velocities, and/or dynamical heating effects), primordial halos can form substantially more numerous and more massive stars. Some of these stars will in turn be accreting more rapidly than they can thermally relax at any given time. Using high-resolution cosmological zoom-in simulations of massive star formation in high-z halos, we find that such rapidly accreting stars produce prominent spectral features that would be detectable by the James Webb Space Telescope. The rapid-accretion episodes within the halo lead to stochastic reprocessing of 0%–20% of the total stellar emission into the rest-frame optical over long timescales, a unique signature which may allow deep observations to identify such objects out to z ∼ 10–13 using mid- and wide-band NIRCam colors alone.

Determining the full satellite population of a Milky Way-mass halo in a highly resolved cosmological hydrodynamic simulation

ABSTRACT We investigate the formation of the satellite galaxy population of a Milky Way-mass halo in a very highly resolved magnetohydrodynamic cosmological zoom-in simulation (baryonic mass resolution mb = 800 $\rm M_{\odot }$). We show that the properties of the central star-forming galaxy, such as the radial stellar surface density profile and star formation history, are (i) robust to stochastic variations associated with the so-called Butterfly Effect and (ii) well converged over 3.5 orders of magnitude in mass resolution. We find that there are approximately five times as many satellite galaxies at this high resolution compared to a standard ($m_b\sim 10^{4-5}\, \rm M_{\odot }$) resolution simulation of the same system. This is primarily because two-thirds of the high-resolution satellites do not form at standard resolution. A smaller fraction (one-sixth) of the satellites present at high-resolution form and disrupt at standard resolution; these objects are preferentially low-mass satellites on intermediate- to low-eccentricity orbits with impact parameters ≲30 kpc. As a result, the radial distribution of satellites becomes substantially more centrally concentrated at higher resolution, in better agreement with recent observations of satellites around Milky Way-mass haloes. Finally, we show that our galaxy formation model successfully forms ultra-faint galaxies and reproduces the stellar velocity dispersion, half-light radii, and V-band luminosities of observed Milky Way and Local Group dwarf galaxies across six orders of magnitude in luminosity (103–$10^{9}\, \rm L_{\odot }$).

Neutral CGM as damped Ly α absorbers at high redshift

ABSTRACT Recent searches for the hosts of z ∼ 4 damped Ly α absorbers (DLAs) have detected bright galaxies at distances of tens of kpc from the DLA. Using the FIRE-2 cosmological zoom simulations, we argue that these relatively large distances are due to a predominantly cool and neutral inner circumgalactic medium (CGM) surrounding high-redshift galaxies. The inner CGM is cool because of the short cooling time of hot gas in ${\lesssim}10^{12}\, {\rm M_{\odot }}$ haloes, which implies that accretion and feedback energy are radiated quickly, while it is neutral due to high volume densities and column densities at high redshift that shield cool gas from photoionization. Our analysis predicts large DLA covering factors (${\gtrsim}50{{\ \rm per\ cent}}$) out to impact parameters ∼0.3[(1 + z)/5]3/2Rvir from the central galaxies at z ≳ 1, equivalent to a proper distance of ${\sim}21\, M_{12}^{1/3} \left(\left(1+z\right)/5\right)^{1/2}\, {\rm kpc}$ (Rvir and M12 are the halo virial radius and mass in units of $10^{12}\, {\rm M_{\odot }}$, respectively). This implies that DLA covering factors at z ∼ 4 may be comparable to unity out to a distance ∼10 times larger than stellar half-mass radii. A predominantly neutral inner CGM in the early universe suggests that its mass and metallicity can be directly constrained by absorption surveys, without resorting to the large ionization corrections as required for ionized CGM.

The interplay of internal and external processes in the buildup of disk galaxies: thick-disk star formation histories in AURIGA simulations

AbstractRecent integral-field spectroscopy observations have revealed that thick- and thin-disk star-formation histories are regulated by the interplay of internal and external processes. We analyze stellar-population properties of 24 spiral galaxies from the AURIGA zoom-in cosmological simulations, to offer a more in-depth interpretation of observable properties. We present edge-on maps of stellar age, metallicity and [Mg/Fe] abundance, and we extract the star-formation and chemical-evolution histories of thin and thick disks. Both show signs of the interplay between internal chemical enrichment and gas and star accretion. Thick disks show particularly complex stellar populations, including an in-situ component, formed from both slowly enriched and accreted more pristine gas, and an additional significant fraction of ex-situ stars.

Satellites around Milky Way Analogs: Tension in the Number and Fraction of Quiescent Satellites Seen in Observations versus Simulations

Abstract We compare the star-forming properties of satellites around Milky Way (MW) analogs from the Stage II release of the Satellites Around Galactic Analogs Survey (SAGA-ii) to those from the APOSTLE and Auriga cosmological zoom-in simulation suites. We use archival GALEX UV imaging as a star formation indicator for the SAGA-ii sample and derive star formation rates (SFRs) to compare with those from APOSTLE and Auriga. We compare our detection rates from the NUV and FUV bands to the SAGA-ii Hα detections and find that they are broadly consistent with over 85% of observed satellites detected in all three tracers. We apply the same spatial selection criteria used around SAGA-ii hosts to select satellites around the MW-like hosts in APOSTLE and Auriga. We find very good overall agreement in the derived SFRs for the star-forming satellites as well as the number of star-forming satellites per host in observed and simulated samples. However, the number and fraction of quenched satellites in the SAGA-ii sample are significantly lower than those in APOSTLE and Auriga below a stellar mass of M * ∼ 108 M ⊙, even when the SAGA-ii incompleteness and interloper corrections are included. This discrepancy is robust with respect to the resolution of the simulations and persists when alternative star formation tracers are employed. We posit that this disagreement is not readily explained by vagaries in the observed or simulated samples considered here, suggesting a genuine discrepancy that may inform the physics of satellite populations around MW analogs.

The AGORA High-resolution Galaxy Simulations Comparison Project. III. Cosmological Zoom-in Simulation of a Milky Way–mass Halo

Abstract We present a suite of high-resolution cosmological zoom-in simulations to z = 4 of a 1012 M ⊙ halo at z = 0, obtained using seven contemporary astrophysical simulation codes (Art-I, Enzo, Ramses, Changa, Gadget-3, Gear, and Gizmo) widely used in the numerical galaxy formation community. The physics prescriptions for gas cooling and heating and star formation are the same as the ones used in our previous Assembling Galaxies of Resolved Anatomy (AGORA) disk comparison but now account for the effects of cosmological processes such as the expansion of the universe, intergalactic gas inflow, and the cosmic ultraviolet background radiation emitted by massive stars and quasars. In this work, we introduce the most careful comparison yet of galaxy formation simulations run by different code groups, together with a series of four calibration steps each of which is designed to reduce the number of tunable simulation parameters adopted in the final run. In the first two steps, we methodically calibrate the gas physics, such as cooling and heating, in simulations without star formation. In the third step, we seek agreement on the total stellar mass produced with the common star formation prescription used in the AGORA disk comparison, in stellar-feedback-free simulations. In the last calibration step, we activate stellar feedback, where each code group is asked to set the feedback prescription to as close to the most widely used one in its code community as possible, while aiming for convergence in the stellar mass at z = 4 to the values predicted by semiempirical models. After all the participating code groups successfully complete the calibration steps, we achieve a suite of cosmological simulations with similar mass assembly histories down to z = 4. With numerical accuracy that resolves the internal structure of a target halo (≲100 physical pc at z = 4), we find that the codes overall agree well with one another, e.g., in gas and stellar properties, but also show differences, e.g., in circumgalactic medium (CGM) properties. We argue that, if adequately tested in accordance with our proposed calibration steps and common parameters, high-resolution cosmological zoom-in simulations can have robust and reproducible results. New code groups are invited to join and enrich this comparison by generating equivalent models or to test the code’s compatibility on their own, by adopting the common initial conditions, the common easy-to-implement physics package, and the proposed calibration steps. Further analyses of the zoom-in simulations presented here will be presented in forthcoming reports from the AGORA Collaboration, including studies of the CGM, simulations by additional codes, and results at lower redshift.