Auriga Project: Milky Way Simulations

• Auriga Project is a suite of high-resolution gravo-magnetohydrodynamical zoom-in simulations designed to study Milky Way-like galaxy formation in a ΛCDM cosmology.
• It employs a moving-mesh code and a zoom-in technique to accurately model baryonic processes, dark matter interactions, and magnetic field dynamics across multiple resolution levels.
• The project provides actionable insights into disc formation, halo structure, circumgalactic gas, and stellar components while offering extensive public data for synthetic survey analyses.

The Auriga Project, introduced by Grand et al., is a suite of cosmological gravo-magnetohydrodynamical “zoom-in” simulations designed to follow the formation and evolution of Milky Way-mass galaxies in a $\Lambda$CDM cosmology with the moving-mesh code AREPO. The original project consists of 30 simulations of isolated Milky Way-mass haloes, while the augmented public-release suite extends this to 40 Milky Way-mass haloes and 26 dwarf galaxy-mass haloes. Across these releases, Auriga has been used to study disc formation, halo structure, circumgalactic gas, stellar haloes, bulges, ultra-diffuse galaxies, magnetic fields, metallicity statistics, and synthetic survey products (Grand et al., 2016, Grand et al., 2024).

1. Scope, cosmology, and numerical design

The original Auriga sample was drawn from a $100\,h^{-1}$ Mpc parent volume, with haloes selected to have final virial masses $1 < M_{200}/10^{12}M_\odot < 2$ and to belong to the most isolated quartile; 30 haloes were then randomly selected from this group. The simulations start at $z=127$ and adopt Planck cosmological parameters, including $\Omega_m = 0.307$, $\Omega_b = 0.048$, $\Omega_\Lambda = 0.693$, and $h = 0.6777$ (Grand et al., 2016).

Auriga employs the zoom-in technique: the high-resolution region is defined by tracing particles within $4\,R_{200}$ of the final halo centre back to the initial conditions, which are sampled from the PANPHASIA Gaussian field. Hydrodynamics and magnetohydrodynamics are solved with AREPO on a dynamic Voronoi mesh, while gravity is computed with a TreePM algorithm. The augmented release describes the code as second-order accurate for self-gravitating MHD in a cosmological context, with Monte Carlo tracer particles available in a subset of runs to follow gas pathways statistically (Grand et al., 2024).

The project is resolved at multiple levels. In the original “level 4” configuration, the baryonic mass resolution is $\sim 5 \times 10^4\,M_\odot$, the dark matter particle mass is $100\,h^{-1}$0, and the force softening reaches a minimum physical value of 369 pc at $100\,h^{-1}$1. The augmented release reports level-4 baryon mass resolution of $100\,h^{-1}$2 and softening of $100\,h^{-1}$3 pc, and level-3 baryon mass resolution of $100\,h^{-1}$4 with softening $100\,h^{-1}$5 pc. The extended suite comprises 40 Milky Way-mass haloes, 12 massive dwarfs, and 14 low-mass dwarfs, with some systems simulated at multiple resolutions (Grand et al., 2016, Grand et al., 2024).

2. Galaxy-formation model and disc-galaxy assembly

Auriga includes a comprehensive galaxy-formation physics model following dark matter, gas, stars, and supermassive black holes. The baryonic model contains primordial and metal-line cooling, a uniform redshift-dependent UV background with self-shielding, a two-phase interstellar medium model, stochastic star formation above a threshold $100\,h^{-1}$6, a star-formation timescale $100\,h^{-1}$7, stellar evolution and chemical enrichment, black-hole growth, AGN feedback, and magnetic fields. Stellar populations adopt a Chabrier IMF (Grand et al., 2016, Grand et al., 2024).

Stellar feedback is implemented with non-local wind particles launched at velocities $100\,h^{-1}$8, where $100\,h^{-1}$9 is the local dark-matter velocity dispersion. Black holes are seeded with $1 < M_{200}/10^{12}M_\odot < 2$0 in haloes with $1 < M_{200}/10^{12}M_\odot < 2$1 and accrete through an Eddington-limited Bondi-Hoyle-Lyttleton prescription supplemented by radio-mode hot-gas accretion; feedback operates in quasar and radio modes. Magnetic fields are seeded at strength $1 < M_{200}/10^{12}M_\odot < 2$2 G comoving at $1 < M_{200}/10^{12}M_\odot < 2$3 and then evolved with ideal MHD (Grand et al., 2016).

A principal result of the original suite is that Auriga reproduces two-component disc-dominated galaxies with appropriate stellar masses, sizes, rotation curves, star formation rates, and metallicities. The simulations also compare well with low-redshift scaling relations including the Tully-Fisher relation, the star-forming main sequence, and HI gas fraction and disc thickness. The star-forming gas discs build rotation and velocity dispersion rapidly for $1 < M_{200}/10^{12}M_\odot < 2$4 before settling into ever-increasing $1 < M_{200}/10^{12}M_\odot < 2$5 ratios (Grand et al., 2016, Grand et al., 2024).

Auriga was specifically designed to investigate disc-size formation mechanisms. Disc scale lengths span 2.2 kpc to 11.6 kpc, and most disc-to-total ratios lie in the range $1 < M_{200}/10^{12}M_\odot < 2$6–$1 < M_{200}/10^{12}M_\odot < 2$7. The dominant driver of present-day disc size is the angular momentum of halo material, with the largest discs produced by quiescent mergers that inspiral into the galaxy and deposit high-angular-momentum material into the pre-existing disc. More violent mergers and strong AGN feedback limit disc size by destroying pre-existing discs and suppressing gas accretion onto the outer disc; the most important factor leading to compact discs is low halo angular momentum (Grand et al., 2016).

Neutral hydrogen provides a complementary view of disc structure. In the vast majority of Auriga systems the HI gas forms extended discs, although more disturbed morphologies are present. The simulations show good agreement with observed HI radial profiles and the HI mass–diameter relation, but the galaxies are systematically larger and more gas-rich than typical nearby galaxies. The amount of HI gas outside the disc plane correlates with star formation rate, consistent with a fountain-like flow interpretation (Marinacci et al., 2016).

3. Dark matter haloes, substructure, and the circumgalactic medium

Auriga has been used extensively to quantify the response of dark-matter structure to baryonic physics. Halo shapes are measured from the reduced inertia tensor,

$1 < M_{200}/10^{12}M_\odot < 2$8

with axis ratios $1 < M_{200}/10^{12}M_\odot < 2$9 and triaxiality parameter

$z=127$0

At $z=127$1, baryons make Milky Way-sized dark-matter haloes rounder at all radii relative to dark-matter-only counterparts. In dark-matter-only runs, haloes are generally prolate in their inner regions with $z=127$2; in full MHD runs they become oblate or nearly spherical, with $z=127$3, especially near the galaxy centre. At distances $z=127$4 kpc, rounder haloes correlate with extended massive stellar discs and low core gas densities. Most galaxies show strong alignment between the halo minor axis and disc angular momentum at all radii, but about 6 of 30 show a substantial radial change of alignment, interpreted as halo twisting. In comparison with stream constraints, about $z=127$5 of the MHD haloes are consistent with the nearly spherical Pal 5 result of Bovy et al. (2016), whereas none of the simulations match the oblate or triaxial Sagittarius-stream constraint of Law and Majewski (2010) (Prada et al., 2019).

Baryonic effects also reshape the subhalo population. Using orbit integration in a time-interpolated axisymmetric potential with GALPY, rather than cubic-spline interpolation, Auriga hydro runs show an 80% reduction in subhalo abundance near the Galactic centre and a 40% reduction within $z=127$6 relative to dark-matter-only simulations for subhaloes with masses $z=127$7–$z=127$8. The corresponding APOSTLE reductions are 50% and 10%, and the difference is traced to the fact that Auriga central galaxies are typically twice as massive as those in APOSTLE. The baryonic suppression persists to $z=127$9 kpc, partly through the destruction of splashback subhaloes (Richings et al., 2018).

The circumgalactic medium in Auriga is likewise highly structured. For 28 $\Omega_m = 0.307$0 Milky Way-mass galaxies, column densities of commonly observed species span $\Omega_m = 0.307$1–$\Omega_m = 0.307$2 dex and covering fractions range from $\Omega_m = 0.307$3 to 90 per cent. The covering fractions of hydrogen and metals correlate positively with stellar mass; the covering fractions of \ion{H}{i}, \ion{C}{iv}, and \ion{Si}{ii} anticorrelate with AGN luminosity because of ionization effects; and the covering fractions of \ion{H}{i}, \ion{C}{iv}, and \ion{Si}{ii} correlate positively with disc fraction, consistent with outflows populating the CGM with cool and dense gas. The study emphasizes that long-term merger assembly history and recent star formation are not the dominant sculptors of the CGM in this sample (Hani et al., 2019).

A direct comparison between matched Auriga and EAGLE zoom simulations isolates the effect of subgrid feedback prescriptions. In this comparison, 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. The difference is attributed to distinct supernova energy-injection methods: in Auriga, gas accretion is almost unaffected by feedback and recycling times are short, whereas EAGLE produces halo-wide outflows that both eject baryons and impede cosmic gas accretion. The paper argues that quasar absorption lines and fast radio burst dispersion measures could discriminate between these regimes (Kelly et al., 2021).

4. Stellar components, accretion histories, and satellite systems

Auriga bulges display a broad range of shapes, sizes, and formation histories, but by standard observational criteria the majority are pseudo-bulges. Photometrically, all Auriga bulges have Sérsic index $\Omega_m = 0.307$4; kinematically, they are rotation-dominated; and none can be classified as a classical bulge. Bulge growth is mostly in situ: the median accreted bulge fraction is $\Omega_m = 0.307$5, 21% of bulges have negligible accreted fraction ($\Omega_m = 0.307$6), and about 75% form most of their in-situ stars inside the present-day bulge region. In 90% of cases, the accreted bulge component originates from fewer than four satellites (Gargiulo et al., 2019).

The stellar haloes of the 30 Milky Way-mass systems are also diverse. Accreted stellar-halo masses span $\Omega_m = 0.307$7 to $\Omega_m = 0.307$8, while the number of significant progenitors contributing 90% of the accreted halo mass ranges from 1 to 14, with a median of 6.5. Halo metallicity gradients, density-profile slopes, ages, and shapes vary substantially, and a key correlation emerges: galaxies with few significant progenitors have more massive haloes, larger negative metallicity gradients, and steeper density profiles. Auriga reproduces the observed stellar-halo mass–metallicity relation, although the simulated haloes are typically too massive (Monachesi et al., 2018).

At lower surface brightness, Auriga identifies 92 ultra-diffuse galaxies whose sizes, central surface brightnesses, Sérsic indices, colours, spatial distribution, and abundance match a wide range of observables. These UDGs inhabit dwarf-mass haloes. Field UDGs arise in low-mass dark-matter haloes with higher spin parameters than normal field dwarfs, while satellite UDGs follow two channels: about 55% formed as field UDGs before accretion, and about 45% were normal field dwarfs that became UDGs through tidal interactions after infall (Liao et al., 2019).

Auriga has also been used to test whether the Local Group satellite mass–metallicity relation excludes severe tidal stripping. For accreted systems, the intrinsic relation between total stellar mass and metallicity is tight, with scatter $\Omega_m = 0.307$9 dex. When only the present-day bound progenitor is considered, the tidally evolved relation has scatter $\Omega_b = 0.048$0 dex, comparable to the observed $\Omega_b = 0.048$1 dex for Milky Way and M31 satellites. Satellites above the relation have typically experienced substantial mass loss; even systems lying exactly on the evolved relation can lose over half of their stellar mass; and only systems substantially below the evolved relation are reliably intact (Riley et al., 8 Sep 2025).

A notable limitation appears in globular-cluster demographics. When all star particles older than 10 Gyr are treated as globular-cluster candidates, Auriga produces candidates at metallicities and radii broadly relevant to the Milky Way system, but overall they are too metal-rich, too radially extended, and do not reproduce the Milky Way’s metallicity bimodality. The run-to-run scatter is smaller than the observed difference between the Milky Way and M31, and the model is judged unlikely to give rise to a globular-cluster system consistent with both galaxies (Halbesma et al., 2019).

5. Magnetic fields and spatial metallicity statistics

Auriga was among the first Milky Way-like cosmological suites to treat magnetic fields self-consistently in a large sample. Starting from a uniform comoving seed field of $\Omega_b = 0.048$2 at $\Omega_b = 0.048$3, the magnetic field is exponentially amplified in halo centres by a small-scale dynamo with an e-folding time of roughly $\Omega_b = 0.048$4. This phase saturates around $\Omega_b = 0.048$5–$\Omega_b = 0.048$6, when the magnetic energy density reaches about $\Omega_b = 0.048$7 of the turbulent energy density and the central field strength is typically $\Omega_b = 0.048$8–$\Omega_b = 0.048$9. At larger radii, differential rotation drives a linear amplification phase that usually saturates between $\Omega_\Lambda = 0.693$0 and $\Omega_\Lambda = 0.693$1. The final radial and vertical field profiles are well described by double exponentials and are in good agreement with observational constraints. The global dynamical impact is small because equipartition is reached only after most star formation has occurred, and the results are numerically well converged (Pakmor et al., 2017).

Chemical inhomogeneity has been studied through two-point metallicity correlation functions in 28 Auriga galaxies. After subtracting the azimuthally averaged metallicity gradient, the correlation function is defined as

$\Omega_\Lambda = 0.693$2

Mock observations of the $\Omega_\Lambda = 0.693$3 snapshots reproduce the correlation functions measured in local galaxy surveys, and the recovered metallicity correlation length $\Omega_\Lambda = 0.693$4 typically matches the true simulation value within $\Omega_\Lambda = 0.693$5. Over cosmic time, individual galaxies show no significant secular evolution in their metallicity correlation functions from $\Omega_\Lambda = 0.693$6 to the present. Fluctuations in $\Omega_\Lambda = 0.693$7 are correlated with, but tend to precede, fluctuations in star formation rate, which the authors interpret as evidence that metal rearrangement occurs at a higher cadence than star formation activity and is more sensitive to mergers, gas inflows and outflows, and fly-bys (Li et al., 2024).

6. Public data products, synthetic observations, and methodological extensions

The augmented Auriga public release provides raw snapshots, group catalogues, merger trees, initial conditions, and supplementary data, together with public analysis tools and worked examples. The released datasets are in HDF5 format and include 128–251 snapshots per run, FOF and SUBFIND catalogues, merger histories, mock catalogues, and high-level tables. Access is provided through Globus, and the accompanying Python package supports snapshot loading, centring and rotating galaxies, reading merger trees, and selecting accreted stars (Grand et al., 2024).

Auriga has been forward-modelled into survey-like stellar catalogues. “Aurigaia” provides mock Gaia DR2 catalogues for six high-resolution Auriga galaxies, with stars selected at $\Omega_\Lambda = 0.693$8 mag and at $\Omega_\Lambda = 0.693$9 mag for $h = 0.6777$0. Two catalogue-generation methods are used: SNAPDRAGONS and a phase-space-preserving method based on Lowing et al. The resulting catalogues include 5-parameter astrometry, radial velocities, multi-band photometry, stellar parameters, dust extinction, uncertainties, and origin flags. Demonstration analyses show that a flared young outer stellar disc should be detectable in A- and B-dwarf stars to radii of $h = 0.6777$1 kpc and that the spin of the stellar halo out to 100 kpc can be accurately measured with Gaia DR2 RR Lyrae stars (Grand et al., 2018).

Auriga has also been processed with SKIRT into dust-aware synthetic images and fluxes from the ultraviolet to the sub-millimetre. Kapoor et al. produced multi-observer datasets for all 30 $h = 0.6777$2 galaxies in 50 broadband filters, with star-forming regions, stellar sources, and diffuse dust treated explicitly. Two dust-assignment recipes were explored: one in which only star-forming or cold gas cells contain dust, and one that assigns dust more broadly to rotationally supported disc gas. Both yield good agreement with observed global dust properties, although Auriga galaxies have slightly higher specific dust masses than observed counterparts. At optical wavelengths the simulated morphologies agree well with observed spiral galaxies, while in the mid- and far-infrared the Auriga galaxies are smaller and more centrally concentrated than the observations (Kapoor et al., 2021).

A more recent methodological development is the “Superstars” scheme, demonstrated for an Auriga Milky Way-like galaxy. This method increases stellar mass resolution by factors of 8 and 64 without changing gas or dark-matter resolution and without altering global properties of the central galaxy or its satellites. It improves the sampling of the stellar disc and halo, reduces numerical heating of the outer disc, keeps substructures in the stellar disc and inner halo more coherent, and makes lower-mass and lower-surface-brightness halo structures more visible. The method is intended for cosmological galaxy simulations that do not resolve individual stars (Pakmor et al., 29 Jul 2025).

Taken together, these results define Auriga as both a simulation suite and a research infrastructure: a numerically converged set of Milky Way-mass and dwarf zoom simulations, a benchmark for baryonic effects on halo structure and the baryon cycle, and a public platform for mock observations and dynamical inference across the disc, halo, and circumgalactic medium (Grand et al., 2016, Grand et al., 2024).