Auriga MHD Cosmological Simulations

• Auriga simulations are high-resolution zoom-in studies that integrate dark matter, gas, stars, black holes, and magnetic fields in a ΛCDM framework.
• They employ the AREPO moving-mesh code with adaptive resolution and sophisticated subgrid physics to capture disc formation, feedback, and merger histories.
• The public data release and detailed analyses provide actionable insights into internal galaxy structure, circumgalactic medium, and baryon cycling.

Auriga is a suite of high-resolution cosmological gravo-magnetohydrodynamical zoom-in simulations designed to follow the coupled evolution of dark matter, gas, stars, supermassive black holes, and magnetic fields in a $\Lambda$CDM universe with the moving-mesh code AREPO. Introduced as a sample of isolated Milky Way-mass haloes and later expanded into an augmented suite including dwarf-galaxy haloes, Auriga has become a common numerical framework for studies of disc formation, bars and boxy/peanut bulges, stellar haloes, dwarf satellites, ultra-diffuse galaxies, the circumgalactic medium, and observationally oriented mock data products (Grand et al., 2016, Grand et al., 2024).

1. Project formation and simulation scope

The original Auriga Project introduced thirty cosmological magneto-hydrodynamical zoom simulations of isolated Milky Way-mass dark haloes. These target haloes were selected from the dark-matter-only counterpart of the EAGLE Ref-L100N1504 cosmological volume, with $1 < M_{200}/(10^{12}\,M_\odot) < 2$, and the final thirty were randomly chosen from the most isolated quartile of candidate haloes. The simulations begin at redshift $127$ (Grand et al., 2016).

The augmented Auriga suite broadened this design beyond the original Milky Way-mass sample. It includes 40 Milky Way-mass haloes with $0.5 < M_{200}/(10^{12}\,M_\odot) < 2$ at $z=0$, 12 dwarf-galaxy-mass haloes with $0.5 < M_{200}/(10^{11}\,M_\odot) < 5$ simulated at both level 4 and level 3, and 14 lower-mass dwarf haloes with $0.5 < M_{200}/(10^{10}\,M_\odot) < 5$ simulated at level 3 only. The initial zoom regions are selected from the EAGLE dark-matter-only simulation, the haloes are chosen from the most isolated quartile to limit environmental contamination, the initial conditions are generated from the Panphasia white-noise field, and the simulations begin at redshift $127$ (Grand et al., 2024).

This progression from a homogeneous Milky Way-mass sample to an expanded Milky Way-plus-dwarf suite is central to the project’s scientific identity. It preserved the original focus on isolated disc-galaxy formation while extending the accessible mass range downward into the dwarf regime and upward into community-scale data release.

2. Numerical architecture and physical model

Auriga is run with AREPO, described in the project literature as a second-order accurate, moving-mesh gravo-magnetohydrodynamics code. AREPO uses TreePM gravity, a Voronoi moving mesh, adaptive resolution that follows the fluid, and finite-volume MHD hydrodynamics. In the H I analysis, the code is further described as quasi-Lagrangian and Galilean-invariant, with second-order Runge–Kutta time integration, least-squares gradient estimates, and Powell’s 8-wave cleaning to control $\nabla\cdot\mathbf{B}=0$ (Grand et al., 2024, Marinacci et al., 2016).

The cosmology is explicitly a $\Lambda$CDM model with Planck 2013 parameters,

$1 < M_{200}/(10^{12}\,M_\odot) < 2$0

with equivalent $1 < M_{200}/(10^{12}\,M_\odot) < 2$1 notation used in several Auriga studies (Grand et al., 2024). Magnetic fields are seeded at the starting redshift with a uniform comoving field

$1 < M_{200}/(10^{12}\,M_\odot) < 2$2

corresponding to a physical strength of $1 < M_{200}/(10^{12}\,M_\odot) < 2$3; the original project paper notes that magnetic fields were a major addition relative to earlier Aquarius-based simulations, although they have little effect on the global stellar disc properties emphasized there (Grand et al., 2024, Grand et al., 2016).

The galaxy-formation model includes primordial and metal-line radiative cooling, UV background heating with self-shielding corrections, a two-phase ISM model for unresolved supernova-regulated gas, stochastic star formation above threshold density $1 < M_{200}/(10^{12}\,M_\odot) < 2$4, stellar evolution, mass loss, and chemical enrichment from AGB stars and SNe Ia/II, stellar feedback with thermal plus kinetic wind energy, supermassive black hole seeding and growth, and AGN feedback in quasar and radio modes (Grand et al., 2024). In the original formulation, the ISM follows the Springel & Hernquist two-phase subgrid model, the star-formation timescale is $1 < M_{200}/(10^{12}\,M_\odot) < 2$5, winds are launched isotropically with

$1 < M_{200}/(10^{12}\,M_\odot) < 2$6

and the wind energy is split equally between kinetic and thermal components; AGN feedback is emphasized as especially important for regulating central star formation and preventing overly massive bulges even in Milky Way-mass systems (Grand et al., 2016).

The project literature repeatedly notes that the Auriga model is broadly similar to IllustrisTNG but differs in the ISM equation of state, AGN feedback implementation, stellar wind treatment, and yield tables (Grand et al., 2024). This combination of a moving-mesh MHD solver and a relatively comprehensive subgrid model is the basis for Auriga’s role as a disc-galaxy formation framework rather than a purely dynamical halo simulation.

3. Resolution strategy and public data release

Auriga uses a small set of well-defined resolution levels. In the augmented suite, level 4 corresponds to baryonic mass resolution $1 < M_{200}/(10^{12}\,M_\odot) < 2$7 and softening $1 < M_{200}/(10^{12}\,M_\odot) < 2$8 pc, while level 3 corresponds to baryonic mass resolution $1 < M_{200}/(10^{12}\,M_\odot) < 2$9 and softening $127$0 pc (Grand et al., 2024). In the original Milky Way-mass runs, the typical high-resolution masses are $127$1 for dark matter and $127$2 for baryons, with the gravitational softening for collisionless particles reaching a maximum physical value of about 369 pc by $127$3 (Grand et al., 2016). These values recur across later analyses of bars, stellar haloes, H I discs, dwarf galaxies, and CGM structure.

A major transition in the project was the first public data release of the augmented suite. That release includes raw snapshots containing gas cells, dark matter particles, stellar/wind particles, and black holes at each output time; group catalogues containing halo and galaxy properties from FoF and SUBFIND; merger trees; full initial-condition files at $127$4, usable with AREPO or GADGET4; supplementary data products including accreted star-particle lists and high-level derived tables; and public analysis tools with worked examples for loading snapshots, reading group catalogues, navigating merger trees, centering and rotating galaxies, and selecting accreted particles (Grand et al., 2024).

The supplementary products extend the release beyond the standard simulation outputs. They include mock Gaia DR2 catalogues for selected simulations, mock Pandas survey catalogues, synthetic UV–submm images and spectral energy distributions from SKIRT, and tabulated high-level data such as stellar disc and bulge properties, H I disc properties, satellite and subhalo statistics, and bar properties. The data are distributed through Globus and the project website, and the analysis code is publicly available via Bitbucket (Grand et al., 2024).

This release materially changed the status of Auriga within numerical galaxy formation. A plausible implication is that the project shifted from a relatively closed simulation program to a reusable reference suite with standardized products for methodological benchmarking, mock-observational pipelines, and direct comparison between independently developed analyses.

4. Disc assembly, internal structure, and gas kinematics

The original Auriga paper demonstrated that the simulations reproduce a wide range of present-day observables, in particular two-component disc-dominated galaxies with appropriate stellar masses, sizes, rotation curves, star-formation rates, and metallicities. The exponential disc scale lengths span

$127$5

and the most important physical conclusion of that work is that disc scale length correlates with halo angular momentum or spin. The largest discs are produced by quiescent mergers that inspiral into the galaxy and deposit high-angular-momentum material into the pre-existing disc, whereas violent major mergers and strong AGN feedback act against disc growth by destroying pre-existing discs or suppressing gas accretion onto the outer disc (Grand et al., 2016).

The augmented suite extended these conclusions into a broader low-redshift and high-redshift context. At $127$6, the simulated galaxies broadly agree with empirical stellar mass–halo mass relations, follow the observed Tully–Fisher relation down to about $127$7, lie on the observed star-forming main sequence from SDSS MPA-JHU, and reproduce the observed H I-to-stellar-mass relation. The same study tracks the kinematics of star-forming gas discs and finds that they build rotation and velocity dispersion rapidly for $127$8, then settle into ever-increasing $127$9. Milky Way-mass systems reach $0.5 < M_{200}/(10^{12}\,M_\odot) < 2$0 by $0.5 < M_{200}/(10^{12}\,M_\odot) < 2$1, present-day values are roughly $0.5 < M_{200}/(10^{12}\,M_\odot) < 2$2 for the heavier Milky Way analogues and $0.5 < M_{200}/(10^{12}\,M_\odot) < 2$3 for the lighter ones, and $0.5 < M_{200}/(10^{12}\,M_\odot) < 2$4 rises from $0.5 < M_{200}/(10^{12}\,M_\odot) < 2$5 at early times to $0.5 < M_{200}/(10^{12}\,M_\odot) < 2$6 at $0.5 < M_{200}/(10^{12}\,M_\odot) < 2$7 for heavier systems and $0.5 < M_{200}/(10^{12}\,M_\odot) < 2$8 for lighter systems, with higher-resolution level 3 Milky Way runs reaching $0.5 < M_{200}/(10^{12}\,M_\odot) < 2$9 at $z=0$0. The absence of quenched massive galaxies in that sample is explicitly noted, possibly linked to AGN mode details and halo selection (Grand et al., 2024).

Auriga has also been used to model the neutral-gas structure of Milky Way-mass discs. Most systems form extended H I discs, often more extended than the stellar discs, with generally good agreement with observed H I radial profiles and the H I mass–diameter relation. At the same time, the Auriga galaxies are systematically larger and more gas-rich than typical nearby galaxies. The amount of H I gas outside the disc plane correlates with the star-formation rate, which is interpreted as consistent with a fountain-like flow (Marinacci et al., 2016).

The suite has become especially prominent in work on non-axisymmetric disc structure. In a cosmological analysis of bar formation and evolution, the $z=0$1 barred fraction is about 60%, declining with redshift to a plateau of $z=0$2 at $z=0$3. Bars are robust and long-lived structures; low- and intermediate-redshift bars form short and then grow, whereas high-redshift bars are born “saturated” in length, likely because of a merger-induced formation pathway. Barred galaxies tend to have lower Toomre $z=0$4 values at the time of formation, are more baryon-dominated at all redshifts, and assemble their stellar mass earlier by about 3 Gyr on average. In the stellar Tully–Fisher relation, barred galaxies are offset from unbarred galaxies by $z=0$5 (Fragkoudi et al., 2024).

Boxy/peanut bulges provide a complementary internal-structure diagnostic. In the dedicated Auriga b/p study, 18 of 30 galaxies are barred at $z=0$6, and 9 of those 18 host a b/p bulge at $z=0$7. B/p bulges form after bars, with a typical delay of $z=0$8–$z=0$9, and every b/p bulge forms after a buckling episode. More than half of the b/p-hosting galaxies show double buckling, Au18 exhibits four buckling episodes, and some systems show temporary weakening or disappearance of the b/p signal before re-emergence. The median size ratio is

$0.5 < M_{200}/(10^{11}\,M_\odot) < 5$0

with a $0.5 < M_{200}/(10^{11}\,M_\odot) < 5$1 median of about 0.48, so the b/p typically extends to about half the bar length (López et al., 19 May 2025).

Bulge structure in Auriga is dominated by pseudo-bulges rather than classical bulges. In the 30-galaxy Milky Way-mass sample, all Auriga bulges have $0.5 < M_{200}/(10^{11}\,M_\odot) < 5$2; the majority are classified as pseudo-bulges, some are composite bulges with a classical component, and none can be classified as a pure classical bulge. Their median accreted fraction is $0.5 < M_{200}/(10^{11}\,M_\odot) < 5$3, 21% have negligible accreted fractions $0.5 < M_{200}/(10^{11}\,M_\odot) < 5$4, and in 90% of bulges the accreted component originates from fewer than four satellites (Gargiulo et al., 2019). A separate observationally motivated analysis of 21 visually barred $0.5 < M_{200}/(10^{11}\,M_\odot) < 5$5 galaxies found realistic bar sizes and luminosities, low-Sérsic bulges, Type I, Type II, and Type III disc breaks, and boxy/peanut structures in 4–6 bars, with Au17 and Au26 identified as buckling systems (Blázquez-Calero et al., 2019).

Auriga has also been used to examine disc stellar kinematics beyond morphology alone. Across 26 disc galaxies, the global stellar-velocity-ellipsoid ratio is

$0.5 < M_{200}/(10^{11}\,M_\odot) < 5$6

and the main conclusion is that similar global SVE values can conceal very different local 2D structures. Bars, fly-bys, and low-mass satellite accretion can imprint strong local signatures even when the global ratio is ordinary (Walo-Martín et al., 2021). In the inner barred regions specifically, $0.5 < M_{200}/(10^{11}\,M_\odot) < 5$7 and $0.5 < M_{200}/(10^{11}\,M_\odot) < 5$8 maps show non-axisymmetric structures that trace the bar light distribution, while $0.5 < M_{200}/(10^{11}\,M_\odot) < 5$9 remains more axisymmetric; the most effective summary statistic is the maximum isophotal difference in the dispersion field (Walo-Martín et al., 2022).

5. Stellar haloes, dwarfs, and the circumgalactic environment

Auriga’s stellar haloes are diverse in their masses, density profiles, metallicities, ages, and shapes, and this diversity is explicitly connected to stochastic accretion and merger histories. In the 30-halo Milky Way-mass suite, the accreted halo mass varies by roughly an order of magnitude, from $0.5 < M_{200}/(10^{10}\,M_\odot) < 5$0 to $0.5 < M_{200}/(10^{10}\,M_\odot) < 5$1; the number of significant progenitors contributing 90% of the accreted halo mass ranges from 1 to 14, with a median of 6.5; and systems with fewer significant progenitors tend to have more massive haloes, larger negative metallicity gradients, and steeper density profiles. The observed halo mass–metallicity relation is reproduced and is set primarily by the dominant satellite contributors (Monachesi et al., 2018).

A notable methodological result from Auriga is that halo metallicity gradients depend strongly on measurement geometry. In a four-galaxy study of Milky Way-like haloes, spherical median $0.5 < M_{200}/(10^{10}\,M_\odot) < 5$2 profiles show strong negative gradients within 100 kpc, but projected minor-axis profiles are significantly flatter. Within the inner 50 kpc, the spherical profile can exceed the minor-axis profile by as much as 0.4 dex, and similar behavior remains even for accreted-only halo stars. This demonstrates that discrepancies between simulated and observed halo metallicity gradients can arise from inconsistent measurement geometry rather than from a fundamental failure of the models (Monachesi et al., 2015).

The dwarf-galaxy applications of Auriga emphasize environmental effects at high redshift. In a filament-based analysis of Auriga dwarfs, galaxies residing in filaments have higher baryonic and stellar fractions than field dwarfs at fixed halo mass, and slightly higher star-formation rates, while no clear overall difference is found in $0.5 < M_{200}/(10^{10}\,M_\odot) < 5$3 colour. At $0.5 < M_{200}/(10^{10}\,M_\odot) < 5$4, the median baryonic fraction in filaments is about five times that in the field, and the gas spin of low-mass dwarf galaxies tends to align with the filament direction at high redshift (Zheng et al., 2021). This supports a picture in which dense, massive high-redshift filaments assist gas accretion and enhance star formation in resident dwarf-sized haloes.

Auriga has also been used to study ultra-diffuse galaxies in a Milky Way-centred environment. Across the 30 level-4 simulations, 92 UDGs are identified: 38 satellites inside the host virial radius and 54 field UDGs between $0.5 < M_{200}/(10^{10}\,M_\odot) < 5$5 and $0.5 < M_{200}/(10^{10}\,M_\odot) < 5$6. Field UDGs are linked to a strong size–spin correlation in low-mass dark matter haloes and have higher median halo spin than blue field normal dwarfs, while satellite UDGs split into two channels: about 55% formed as field UDGs before accretion, and about 45% were normal field dwarfs transformed into UDGs by tidal interactions after infall. Their mean halo mass is $0.5 < M_{200}/(10^{10}\,M_\odot) < 5$7, with a maximum halo mass of $0.5 < M_{200}/(10^{10}\,M_\odot) < 5$8, which places them in the dwarf regime rather than the “failed $0.5 < M_{200}/(10^{10}\,M_\odot) < 5$9” regime (Liao et al., 2019).

The circumgalactic medium around Auriga galaxies is likewise highly diverse. For 28 selected $127$0 Milky Way-mass haloes, commonly observed ionic column densities span about 3–4 dex and covering fractions range from $127$1 to 90%. The covering fractions of the parent elements H, C, O, and Si correlate positively with stellar mass; the covering fractions of H I, C IV, and Si II anticorrelate with AGN bolometric luminosity because of ionization effects; and the covering fractions of H I, C IV, and Si II correlate positively with disc fraction, while no significant correlation is found with 10 Myr- or 500 Myr-averaged specific star-formation rate or with ex-situ stellar mass fraction as a proxy for long-term merger history (Hani et al., 2019). The implication drawn in that work is that halo mass and recent star formation alone are insufficient descriptors of Milky Way-mass CGM structure.

6. Comparative interpretation, methodological cautions, and legacy

A central lesson from Auriga is that realistic stellar galaxies do not uniquely determine a unique baryon cycle. In a controlled Local Group-like comparison, the same zoom-in initial conditions were evolved with the Auriga and EAGLE galaxy-formation models. Auriga predicts that the Milky Way is almost baryonically closed, whereas EAGLE predicts that only half of the expected baryons reside within the halo and that the baryon deficiency extends to the Local Group scale. The study attributes the difference primarily to the distinct energy-injection methods from supernova feedback to gas, with Auriga’s feedback behaving in a more fountain-like, recycling-dominated manner (Kelly et al., 2021). This establishes subgrid sensitivity as a substantive interpretive issue rather than a minor implementation detail.

Auriga’s realism is therefore heterogeneous rather than uniform across all diagnostics. In globular-cluster-system tests, the simulations produce enough old stellar mass for GC candidates, but those age-selected candidates are too metal-rich overall and generally found at larger radii than observed in the Milky Way and M31. The model is judged unlikely to give rise to a globular-cluster system consistent with both galaxies simultaneously (Halbesma et al., 2019). In the augmented Auriga overview, the absence of quenched massive galaxies in the sample is explicitly noted as a possible consequence of AGN mode details and halo selection (Grand et al., 2024). These results do not negate the suite’s successes, but they delimit the regimes in which direct observational agreement can be claimed without qualification.

Auriga has also exposed limits in widely used analysis methodologies. A 2025 study testing integrals-of-motion clustering in Solar-neighbourhood-like halo samples found that the method is very effective for debris from accretion events less than about 6–7 Gyr old, but struggles to detect most debris from older accretion. Most detected structures suffer from significant contamination by in-situ stars, and the background-halo construction can generate artificial detections (Thomas et al., 14 Apr 2025). When read alongside the metallicity-geometry result for stellar haloes, this suggests that Auriga’s long-term importance lies not only in producing simulated galaxies, but also in furnishing controlled counterexamples to overly direct observational inference.

Taken together, the Auriga magneto-hydrodynamical simulations define a broad numerical program in which Milky Way-mass and dwarf galaxies are studied across cosmic time with a common MHD zoom-in framework, common cosmology, and common subgrid model, but with enough halo-to-halo diversity to support population-level inference. Their most durable contribution is the combination of detailed internal galaxy structure, cosmological assembly history, and now a public data ecosystem, allowing the same simulation suite to support both physical interpretation and methodological audit across a large fraction of contemporary disc-galaxy research.