HESTIA Local Group Simulations

• HESTIA simulations are high-resolution constrained cosmological models that reproduce the Local Group’s structure by using observed peculiar velocity fields.
• They employ a hierarchical zoom strategy with advanced magneto-hydrodynamic galaxy formation physics to accurately resolve Milky Way and Andromeda analogues and surrounding cosmic features.
• Results validate realistic observables like rotational curves, satellite distributions, and stellar mass functions, providing actionable insights into galaxy evolution and dynamics.

HESTIA simulations are a suite of constrained cosmological simulations of the Local Group, designed to form a Milky Way–Andromeda system inside a realistically reconstructed nearby cosmic environment. The acronym is given as “High-resolutions Environmental Simulations of The Immediate Area.” Their defining feature is the use of initial conditions constrained by the observed peculiar velocity field of nearby galaxies, so that by $z=0$ the simulations reproduce the local cosmography, including the Local Void, the Local Sheet, the Local Filament, and a Virgo-like cluster, while simultaneously resolving the Milky Way and M31 analogues with magneto-hydrodynamic galaxy-formation physics (Libeskind et al., 2020).

1. Constrained Local-Universe construction

HESTIA was developed within the CLUES framework as a Local-Group-specific realization of constrained $\Lambda$CDM structure formation. The parent simulations use Wiener-Filter/Constrained-Realization reconstructions of the CosmicFlows-2 catalog, described in the project paper as a catalog of $\sim 8\,000$ galaxy peculiar velocities, in order to impose the observed large-scale density and velocity field on the initial conditions (Libeskind et al., 2020). In the notation used for the constrained realizations, the initial density field may be written as

$$\delta_{\rm initial} = \delta_{\rm WF} + \delta_{\rm CR},$$

where $\delta_{\rm WF}$ is the Wiener-filtered field and $\delta_{\rm CR}$ is a constrained random residual consistent with the data (Dupuy et al., 2022).

The selection strategy proceeds hierarchically. A low-resolution dark-matter-only ensemble of about $1\,000$ constrained realizations is evolved to $z=0$, and candidate Local Groups are selected by explicit cosmographic cuts. In the formulation given for the original suite, the candidate pair must lie within $5$ Mpc of the origin; each halo must satisfy $$M_{\rm halo}\in[8\times10^{11},3\times10^{12}]\,{\rm M_\odot}$$; the pair separation must be $\Lambda$0 Mpc; no third halo more massive than the smaller primary may lie within $\Lambda$1 Mpc; the mass ratio must be at least $\Lambda$2; and the relative radial velocity must be negative. The same selection also requires a Virgo Cluster analogue with mass at least $\Lambda$3, located within $\Lambda$4 Mpc of the observed supergalactic coordinates, the Virgo–LG distance to agree within $\Lambda$5 Mpc of the true value, and no cluster more massive than Virgo within $\Lambda$6 Mpc (Libeskind et al., 2020).

This construction makes HESTIA distinct from unconstrained zoom simulations of isolated halo pairs. The suite is explicitly intended to preserve both the internal properties of the MW–M31 system and the external tidal field of the Local Volume. A plausible implication is that HESTIA is best understood not merely as a pair simulation, but as a Local Group simulation embedded in a reconstructed cosmographic context.

2. Numerical realization and baryonic physics

The suite uses a multi-tier zoom strategy. Low-resolution candidate selection is followed by intermediate- and high-resolution resimulations of the accepted Local Group environments. The core numerical configurations reported for the project are summarized below (Libeskind et al., 2020).

Configuration	High-resolution region	Resolution
DM-only precursor	central $\Lambda$7 Mpc ($\Lambda$8 Mpc $\Lambda$9) region at $\sim 8\,000$0 effective resolution	$\sim 8\,000$1
Intermediate zoom	$\sim 8\,000$2 Mpc ($\sim 8\,000$3 Mpc $\sim 8\,000$4) sphere around the LG center	$\sim 8\,000$5, $\sim 8\,000$6, $\sim 8\,000$7 pc
High-resolution zoom	two overlapping $\sim 8\,000$8 Mpc ($\sim 8\,000$9 Mpc $$\delta_{\rm initial} = \delta_{\rm WF} + \delta_{\rm CR},$$0) spheres around the two main haloes	$$\delta_{\rm initial} = \delta_{\rm WF} + \delta_{\rm CR},$$1, $$\delta_{\rm initial} = \delta_{\rm WF} + \delta_{\rm CR},$$2, $$\delta_{\rm initial} = \delta_{\rm WF} + \delta_{\rm CR},$$3 pc

The hydrodynamics is solved with the moving-mesh code AREPO on a quasi-Lagrangian Voronoi mesh, using ideal MHD and a hybrid Tree-PM gravity solver (Libeskind et al., 2020). The project paper explicitly gives the ideal-MHD system solved by the code, including continuity,

$$\delta_{\rm initial} = \delta_{\rm WF} + \delta_{\rm CR},$$4

momentum,

$$\delta_{\rm initial} = \delta_{\rm WF} + \delta_{\rm CR},$$5

and induction,

$$\delta_{\rm initial} = \delta_{\rm WF} + \delta_{\rm CR},$$6

The magnetic field evolution uses the 8-wave Powell scheme, and the Riemann solver is HLLD in the frame of each moving face (Libeskind et al., 2020).

Galaxy formation follows the Auriga model. The subgrid modules include primordial and metal-line cooling, a uniform UV background, a two-phase ISM, stochastic star formation above $$\delta_{\rm initial} = \delta_{\rm WF} + \delta_{\rm CR},$$7, stellar evolution, chemical enrichment, supernova-driven winds, SMBH seeding and growth, AGN feedback, and a uniform seed magnetic field $$\delta_{\rm initial} = \delta_{\rm WF} + \delta_{\rm CR},$$8 G at $$\delta_{\rm initial} = \delta_{\rm WF} + \delta_{\rm CR},$$9 (Libeskind et al., 2020). Later HESTIA papers also describe the production runs as using the Auriga-flavored model with radiative cooling and heating, star formation, stellar feedback, black-hole growth, AGN feedback, and ideal MHD (Dupuy et al., 2022).

The cosmology reported for the main suite is Planck-compatible, with $\delta_{\rm WF}$0, $\delta_{\rm WF}$1, $\delta_{\rm WF}$2, $\delta_{\rm WF}$3, and $\delta_{\rm WF}$4 (Libeskind et al., 2020).

3. Local Group analogues and empirical fidelity

The central claim of HESTIA is that the constrained environment does not prevent the formation of realistic Local Group primaries. In the project paper, the simulated MW and M31 analogues have halo masses $\delta_{\rm WF}$5–$\delta_{\rm WF}$6 and $\delta_{\rm WF}$7–$\delta_{\rm WF}$8, a mass ratio $\delta_{\rm WF}$9–$\delta_{\rm CR}$0, separations of $\delta_{\rm CR}$1–$\delta_{\rm CR}$2 kpc, radial velocities from about $\delta_{\rm CR}$3 to $\delta_{\rm CR}$4 km s$\delta_{\rm CR}$5, and tangential velocities of about $\delta_{\rm CR}$6–$\delta_{\rm CR}$7 km s$\delta_{\rm CR}$8 (Libeskind et al., 2020). Three named high-resolution realizations—09_18, 17_11, and 37_11—are repeatedly used in subsequent analyses (Salomon et al., 2023).

Several standard observables were used to validate the primaries. The circular-velocity curves

$\delta_{\rm CR}$9

match the approximately flat Milky Way curve at $1\,000$0 km s$1\,000$1 and the M31 curve at $1\,000$2 km s$1\,000$3. Stellar masses lie on the empirical Guo et al. relation, and the surface-brightness profiles are described by a Sersic bulge plus an exponential disc,

$1\,000$4

with fitted values $1\,000$5–$1\,000$6 kpc, $1\,000$7–$1\,000$8 kpc, $1\,000$9–$z=0$0, and disc-to-total ratios $z=0$1–$z=0$2 (Libeskind et al., 2020).

Satellite populations constitute a second validation axis. In HESTIA, well-resolved satellites within $z=0$3 kpc reproduce the observed stellar mass functions and $z=0$4-band luminosity functions down to $z=0$5 and $z=0$6, and the cumulative radial satellite distributions broadly match those of the Milky Way and M31 out to $z=0$7 kpc (Libeskind et al., 2020). The most massive-satellite statistics also admit Magellanic analogues: about $z=0$8 of hosts have a satellite with $z=0$9, with some located at $5$0 kpc and others at $5$1–$5$2 kpc (Libeskind et al., 2020).

The suite also reports a specific environmental effect on assembly histories. Constrained Local Group haloes assemble half of their $5$3 mass at $5$4, about $5$5 Gyr earlier than unconstrained analogues (Libeskind et al., 2020). This suggests that the local tidal field and filamentary inflow are not merely boundary conditions, but active determinants of the detailed growth history.

4. Cosmic-web feeding, barycentric infall, and multiphase gas

A major scientific use of HESTIA has been the study of how the Local Group is fed by the surrounding cosmic web. In the accretion analysis of the three high-resolution runs, satellites are identified with AHF and tracked back to their first crossing of $5$6, where

$5$7

The results show two accretion eras separated by a trough around $5$8, a median pre-infall travel distance of $5$9 Mpc with a tail to $$M_{\rm halo}\in[8\times10^{11},3\times10^{12}]\,{\rm M_\odot}$$0 Mpc, and a strong alignment between infall directions and the $$M_{\rm halo}\in[8\times10^{11},3\times10^{12}]\,{\rm M_\odot}$$1 axis of both the tidal and velocity-shear tensors. The alignment is strongest for the early-infall population with $$M_{\rm halo}\in[8\times10^{11},3\times10^{12}]\,{\rm M_\odot}$$2 (Dupuy et al., 2022).

HESTIA has also been used to interpret Local Group gas and galaxy kinematics in observational reference frames. From a Sun-like observer placed in the simulated Milky Way disc, the simulated sky maps show that after subtracting galactic rotation, material outside the Milky Way virial radius develops a radial-velocity dipole aligned with the Local Group barycentre direction when the MW–M31 radial velocity is close to the observed value. In the 17_11 run, where the pair has $$M_{\rm halo}\in[8\times10^{11},3\times10^{12}]\,{\rm M_\odot}$$3 km s$$M_{\rm halo}\in[8\times10^{11},3\times10^{12}]\,{\rm M_\odot}$$4, the GSR-frame dipole amplitude is of order $$M_{\rm halo}\in[8\times10^{11},3\times10^{12}]\,{\rm M_\odot}$$5–$$M_{\rm halo}\in[8\times10^{11},3\times10^{12}]\,{\rm M_\odot}$$6 km s$$M_{\rm halo}\in[8\times10^{11},3\times10^{12}]\,{\rm M_\odot}$$7, and transformation to the LGSR narrows the velocity histograms and produces the same qualitative dipole sharpening seen in absorption-line data (Biaus et al., 2022).

The multiphase circumgalactic and intragroup medium is another recurrent HESTIA theme. Mock skymaps of ion columns show that low-temperature tracers such as H I and Si III are more clumpy than O VI, O VII, and O VIII, and the angular power spectra confirm that small-scale power is highest for H I, next for Si III, and lowest for O VI, O VII, and O VIII. HESTIA under-produces the M31 Project AMIGA columns, but remains consistent with low-redshift galaxy samples; one proposed explanation is contamination of M31 sightlines by Milky Way CGM gas (Damle et al., 2022).

A later ion-by-ion kinematic analysis extends this Local Group picture. In two highest-resolution runs, H I, Si III, and C IV primarily trace cold gas inside the Milky Way and Andromeda haloes, whereas O VI, O VII, and O VIII trace hot halo and intragroup gas. After filtering out disc-like rotation, sightlines toward the barycentre are more likely to be dominated by external Local Group material, and the pressure outside the Milky Way halo is systematically higher toward the barycentre direction than toward its antipode (Biaus et al., 19 Jan 2026). A plausible implication is that HESTIA treats the Local Group not as two disconnected CGM reservoirs, but as a dynamically coupled multiphase medium.

5. Stellar-halo assembly and dynamical non-equilibrium

HESTIA has been used extensively for stellar-halo archaeology. In the analysis of the in-situ component, each of the six MW/M31 analogues experiences between one and four mergers with stellar mass ratios between $$M_{\rm halo}\in[8\times10^{11},3\times10^{12}]\,{\rm M_\odot}$$8 and $$M_{\rm halo}\in[8\times10^{11},3\times10^{12}]\,{\rm M_\odot}$$9 relative to the host at the time of the merger, with one exception these significant events occurring $\Lambda$00–$\Lambda$01 Gyr ago. The inner stellar halo contains an in-situ fraction of about $\Lambda$02–$\Lambda$03, while beyond $\Lambda$04 kpc this fraction typically does not exceed $\Lambda$05. Significant mergers sharply increase the orbital eccentricity and reduce the rotational velocity $\Lambda$06 of pre-existing disc stars, reproducing Splash- and Plume-like features, and the $\Lambda$07 plane develops wedge structures mainly populated by stars born between significant mergers (Khoperskov et al., 2022).

The accreted component shows equally strong complexity. Across the same six galaxies, there are a few dozen mergers in total, but only $\Lambda$08–$\Lambda$09 with stellar mass ratio $\Lambda$10; depending on the halo definition, the most massive merger contributes between $\Lambda$11 and $\Lambda$12 of the total stellar halo. Individual merger remnants overlap strongly in $\Lambda$13–$\Lambda$14, $\Lambda$15–$\Lambda$16, and $\Lambda$17 space, and their loci move with time because the host mass grows and the potential is non-axisymmetric. All six galaxies reveal radially hot, non-rotating or weakly counter-rotating Gaia-Sausage-like structures in the $\Lambda$18–$\Lambda$19 plane (Khoperskov et al., 2022).

The chemical-abundance analysis adds a third dimension to this reconstruction. Accreted debris are chemically distinct from surviving dwarf galaxies, and accreted stellar haloes reveal abundance gradients in $\Lambda$20, where the most metal-rich stars formed in the inner parts of disrupted systems and contribute preferentially to the central host regions. Prograde accreted stars exhibit a prominent knee in the $\Lambda$21–$\Lambda$22 plane, whereas retrograde accreted stars typically occupy a high-$\Lambda$23 sequence. At $\Lambda$24, the in-situ metal-poor stars show between zero and $\Lambda$25 km s$\Lambda$26 net rotation, consistent with an Aurora-like population (Khoperskov et al., 2022).

Several later HESTIA studies show that the present-day Local Group analogues are dynamically non-equilibrium systems. In the centre-of-mass analysis, the all-particle COM is dominated by dark matter but closely tracked by stars; in quiescent hosts the velocity offsets are marginal, $\Lambda$27 km s$\Lambda$28, while total COM–disc positional offsets at $\Lambda$29 span about $\Lambda$30–$\Lambda$31 kpc, exceeding $\Lambda$32 kpc in half the cases. In one MW analogue, a recent $\Lambda$33 satellite at $\Lambda$34 kpc drives $\Lambda$35 km s$\Lambda$36 at $\Lambda$37 (Salomon et al., 2023). In an allied analysis of the time-dependent gravitational field, the total potential is expanded as

$\Lambda$38

with $\Lambda$39 and radial splines out to $\Lambda$40 Mpc. Over the last $\Lambda$41 Gyr, the anisotropic coefficients vary by order $\Lambda$42–$\Lambda$43, and restricting the expansion to $\Lambda$44 underestimates the quadrupole at $\Lambda$45 kpc by at least $\Lambda$46 (Arakelyan et al., 2024). This suggests that HESTIA supports a non-stationary description of the Milky Way potential at large radii.

6. Derived predictions and later extensions

Because HESTIA resolves the Local Group in its environmental context, it has been used to make forward predictions for several nearby-galaxy and circumgalactic observables. One example is the predicted population of Local Group ultra-diffuse galaxies. For a Local Group with enclosed mass $\Lambda$47, the hydrodynamic constrained simulations yield a forecast of $\Lambda$48 isolated UDGs within $\Lambda$49 Mpc, for $\Lambda$50 and $\Lambda$51 kpc. Of these, $\Lambda$52 are expected to be detectable in the footprint of SDSS, while an all-sky survey with SDSS-, DES-, or LSST-like depth would observe almost the entire Local Group field population (Newton et al., 2022).

A second extension concerns dwarf-galaxy circumgalactic media. In a Magellanic-analog pair found in HESTIA, the massive dwarf has $\Lambda$53 and hosts a warm coronal envelope with $\Lambda$54 K, while the interacting companion system produces a neutral H I stream extending over $\Lambda$55 kpc with $\Lambda$56. Surveying the suite, all halos with $\Lambda$57 host warm coronae at $\Lambda$58 (Chisholm et al., 21 Apr 2025).

HESTIA has also been used as the Local Group component of dispersion-measure models for fast radio bursts. In one implementation based on the 37_11 high-resolution run, the Milky Way halo contribution excluding the NE2001 disc spans $\Lambda$59–$\Lambda$60 pc cm$\Lambda$61 with mean $\Lambda$62 pc cm$\Lambda$63 and $\Lambda$64 pc cm$\Lambda$65, while the intragroup medium between $\Lambda$66 and $\Lambda$67 Mpc contributes $\Lambda$68 pc cm$\Lambda$69 with $\Lambda$70 (Huang et al., 2024).

Radiative-transfer post-processing has further extended HESTIA into the reionization epoch. In the reionization study based on the 09_18 realization, a uniform $\Lambda$71 dark-matter run is used to calibrate source models, which are then applied to a $\Lambda$72-effective zoom resolving haloes down to $\Lambda$73. In all scenarios, reionization of the Local Group proceeds in an inside-out manner; the MW and M31 progenitors reach $\Lambda$74 ionization at $\Lambda$75–$\Lambda$76, earlier than the global midpoint at $\Lambda$77–$\Lambda$78, and external ionization fronts play a negligible role. Present-day satellite reionization redshifts show only a weak correlation with present-day host distance, while satellites assembled before reionization are systematically more massive today (Attard et al., 12 Sep 2025).

Taken together, these extensions show that HESTIA has evolved from a constrained Local Group formation suite into a multipurpose numerical laboratory for nearby-galaxy structure, Local Group gas dynamics, satellite archaeology, non-equilibrium gravitational modeling, and observable forecasting. A plausible implication is that its principal scientific value lies in combining Local Volume realism with zoom-level baryonic resolution, rather than treating either of those design goals in isolation.