Hestia Constrained Simulations

• Hestia constrained simulations are high-resolution cosmological models that use observed peculiar velocities to recreate a Milky Way–Andromeda analogue within the local cosmic web.
• They employ a combination of Wiener-filter reconstruction, Constrained Realizations, and the Reverse Zel’dovich Approximation to accurately set initial conditions.
• The simulations enable detailed studies of local group kinematics, environmental influences, and galaxy formation processes in a realistic, constrained cosmographic setting.

Hestia constrained simulations are a set of cosmological simulations of the Local Group in which initial conditions constrained by the observed peculiar velocity of nearby galaxies are employed to simulate the local cosmography and to form a Milky Way–Andromeda analogue pair embedded in the surrounding Virgo cluster, Local Void, and Local Filament (Libeskind et al., 2020). Within the HESTIA framework, the constrained large-scale realization is built by combining Wiener-filter reconstruction, constrained realizations, and the Reverse Zel’dovich Approximation (RZA), after which selected Local Group candidates are re-simulated at high resolution with magneto-hydrodynamics and Auriga galaxy-formation physics (Doumler et al., 2012).

1. Constrained-realization basis

The HESTIA initial conditions are rooted in the CLUES framework and use peculiar-velocity data from CosmicFlows-2, typically after grouping to suppress virial motions and applying bias-minimization techniques, so that the reconstructed large-scale density and velocity fields reproduce the observed Local Volume (Dupuy et al., 2022). In the linear regime, the peculiar velocity and density perturbation fields are related through

$$\mathbf{v}(\mathbf{k}) = i\,H_0\,f\,\delta(\mathbf{k})\,\mathbf{k}/k^2,$$

or equivalently

$$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$

with $\nabla^2 \Phi(\mathbf{x})=\delta(\mathbf{x})$ (Libeskind et al., 2020).

The constrained-realization step imposes the observational information while restoring the missing small-scale power. In Fourier space, the mean constrained mode may be written as

$\delta_{\rm CR}(k)=P(k)\,[P(k)+N(k)]^{-1}d(k),$

where $P(k)$ is the prior matter power spectrum, $N(k)$ is the noise power, and $d(k)$ is the Fourier transform of the data vector (Doumler et al., 2012). In real space, one may also express the correction to a random realization $\delta_{\rm RR}$ as

$$\delta_{\rm CR}(\mathbf{x})=\delta_{\rm RR}(\mathbf{x})+\sum_{i,j}\langle\delta(\mathbf{x})\,c_i\rangle\,[\langle c_i\,c_j\rangle]^{-1}(c_j-\tilde c_j),$$

where $\tilde c_j$ is the value of constraint $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$0 measured in $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$1 (Doumler et al., 2012).

A defining feature of the HESTIA-style pipeline is the Lagrangian augmentation introduced by RZA. Linear theory gives the Zel’dovich mapping

$$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$2

so a first-order reverse estimate of the displacement from the Wiener-filtered velocity field is

$$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$3

Each observational constraint is therefore relocated to a pseudo-Lagrangian coordinate

$$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$4

and the original velocity constraints are applied at $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$5 rather than at $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$6 (Doumler et al., 2012). This removes the large, $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$7 cosmic shift that affected earlier constrained-realization-only reconstructions (Doumler et al., 2012).

2. From peculiar velocities to HESTIA initial conditions

The practical workflow starts from radial peculiar velocities $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$8 of galaxies, together with errors $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$9 and angular positions. Galaxies are combined into groups to suppress virial motions and retain the large-scale coherent flow. The noise covariance is then built as $\nabla^2 \Phi(\mathbf{x})=\delta(\mathbf{x})$0 plus a small non-linearity floor to avoid singularities (Doumler et al., 2012).

From the grouped radial data, one computes the Wiener-filtered three-dimensional velocity field,

$\nabla^2 \Phi(\mathbf{x})=\delta(\mathbf{x})$1

forms the RZA shift $\nabla^2 \Phi(\mathbf{x})=\delta(\mathbf{x})$2, and relocates each constraint from $\nabla^2 \Phi(\mathbf{x})=\delta(\mathbf{x})$3 to $\nabla^2 \Phi(\mathbf{x})=\delta(\mathbf{x})$4 while keeping the measured velocity component unchanged (Doumler et al., 2012). The constrained density mode $\nabla^2 \Phi(\mathbf{x})=\delta(\mathbf{x})$5 is then constructed using these shifted constraints.

The final particle initial conditions follow from

$\nabla^2 \Phi(\mathbf{x})=\delta(\mathbf{x})$6

after which dark-matter particles are placed initially on a regular grid or glass and displaced to

$\nabla^2 \Phi(\mathbf{x})=\delta(\mathbf{x})$7

with particle velocities assigned as

$\nabla^2 \Phi(\mathbf{x})=\delta(\mathbf{x})$8

If hydrodynamics is included, each grid cell’s mass is split into dark matter and gas according to $\nabla^2 \Phi(\mathbf{x})=\delta(\mathbf{x})$9, and the same displacement and velocity fields are applied to both components (Doumler et al., 2012).

In HESTIA, this constrained reconstruction is coupled to a screening stage in which an ensemble of low-resolution dark-matter-only runs in a $\delta_{\rm CR}(k)=P(k)\,[P(k)+N(k)]^{-1}d(k),$0 box is searched for suitable Local Group analogues. Typical selection criteria require two halos with $\delta_{\rm CR}(k)=P(k)\,[P(k)+N(k)]^{-1}d(k),$1, separation $\delta_{\rm CR}(k)=P(k)\,[P(k)+N(k)]^{-1}d(k),$2–$\delta_{\rm CR}(k)=P(k)\,[P(k)+N(k)]^{-1}d(k),$3, mass ratio $\delta_{\rm CR}(k)=P(k)\,[P(k)+N(k)]^{-1}d(k),$4, no comparably massive third neighbor within $\delta_{\rm CR}(k)=P(k)\,[P(k)+N(k)]^{-1}d(k),$5, and approaching radial velocity (Biaus et al., 2022). Once such pairs are found, the central region is re-simulated at progressively higher resolution (Biaus et al., 2022).

3. Numerical architecture of the HESTIA suite

The parent HESTIA volume is a constrained $\delta_{\rm CR}(k)=P(k)\,[P(k)+N(k)]^{-1}d(k),$6 box. Initial low-resolution runs use $\delta_{\rm CR}(k)=P(k)\,[P(k)+N(k)]^{-1}d(k),$7 dark-matter particles, while the Local Group region is refined first to an effective $\delta_{\rm CR}(k)=P(k)\,[P(k)+N(k)]^{-1}d(k),$8 resolution and then, in the highest-resolution realizations, to overlapping $\delta_{\rm CR}(k)=P(k)\,[P(k)+N(k)]^{-1}d(k),$9 zoom regions around the two main halos (Libeskind et al., 2020). One formulation describes a $P(k)$0 sphere at $P(k)$1 effective resolution, followed by two overlapping $P(k)$2 spheres around the Milky Way and M31 haloes refined to $P(k)$3 (Biaus et al., 19 Jan 2026).

The high-resolution runs employ the moving-mesh Voronoi code AREPO together with the Auriga galaxy-formation model (Biaus et al., 19 Jan 2026). The highest-resolution zooms have dark-matter particle mass $P(k)$4, gas-cell target mass $P(k)$5, and gravitational softening $P(k)$6 in the high-resolution region (Khoperskov et al., 2022). The cosmology used across the main HESTIA suite is $P(k)$7, $P(k)$8, $P(k)$9, $N(k)$0, and $N(k)$1 (Biaus et al., 2022).

The physical model includes ideal MHD, primordial and metal-line cooling, a UV background, star formation in a multiphase interstellar medium, stellar evolution, chemical enrichment, supernova feedback, black-hole seeding and growth, AGN feedback, and isotropic, mass-loading-driven stellar winds (Biaus et al., 2022). In one detailed description, the seed magnetic field is $N(k)$2 at $N(k)$3 (Biaus et al., 2022). The resulting suite contains three principal high-resolution realizations, identified as 09_18, 17_11, and 37_11, each containing two dominant halos with $N(k)$4–$N(k)$5 (Biaus et al., 2022).

4. Accuracy, reconstruction performance, and Local Group fidelity

The methodological benchmark for HESTIA-style constrained simulations comes from the mock-data tests of the RZA+CR procedure. In those tests, the reconstructed initial conditions were evolved forward to $N(k)$6 and compared against the original reference simulation. The addition of RZA to the constrained-realization method significantly improved both the reconstruction of the initial conditions and the accuracy of the resimulations (Doumler et al., 2012). Haloes from the original simulation were recovered with an average positional accuracy of about $N(k)$7 and a factor of $N(k)$8 in mass, down to haloes with mass $N(k)$9, whereas CR-only reconstructions recovered only the most massive haloes, about $d(k)$0 and higher, with a systematic positional shift of about $d(k)$1 due to the cosmic displacement field (Doumler et al., 2012).

At the level of the full HESTIA suite, the constrained environment is designed to recover the Virgo cluster, Local Void, and Local Filament, together with the mutual approach of the Milky Way and M31 (Biaus et al., 19 Jan 2026). The three principal high-resolution realizations span present-day MW–M31 separations of $d(k)$2–$d(k)$3 and relative radial velocities of $d(k)$4, $d(k)$5, and $d(k)$6 for seeds 37_11, 9_18, and 17_11, respectively, with 17_11 being the closest to the observed $d(k)$7 (Biaus et al., 2022).

The Hestia project reported that the simulated Local Group galaxies resemble the Milky Way and Andromeda in halo mass, mass ratio, stellar disc mass, morphology separation, relative velocity, rotation curves, bulge-disc morphology, satellite galaxy stellar mass function, satellite radial distribution, and in some cases the presence of a Magellanic cloud like object (Libeskind et al., 2020). This suggests that the constrained large-scale environment is not merely a boundary condition, but an active component of the model’s success in reproducing Local Group-scale structure and kinematics.

5. Scientific results enabled by the constrained environment

A substantial part of the HESTIA literature uses the constrained Local Group context to study multiphase gas. The $d(k)$8 gas around the Milky Way–M31 system shows that low-temperature tracers such as H I and Si III are more clumpy than warmer tracers such as O VI, O VII, and O VIII, and that HESTIA under-produces the column densities of the M31 observations while remaining consistent with observations of low-redshift galaxies (Damle et al., 2022). A proposed explanation is contamination of the M31 measurements by gas in the circumgalactic medium of the Milky Way (Damle et al., 2022).

In Local Group kinematics, HESTIA produces a radial-velocity dipole once Galactic rotation is removed. In the approaching realizations 9_18 and 17_11, sightlines toward the barycentre are on average approaching in the Galactic Standard of Rest, while those in the opposite direction are lagging; the effect vanishes in 37_11, where the pair is receding (Biaus et al., 2022). A pure-dipole description gives $d(k)$9–$\delta_{\rm RR}$0, roughly aligned with M31, and in the best-matched run 17_11, $\delta_{\rm RR}$1 with $\delta_{\rm RR}$2 (Biaus et al., 2022). A later ion-resolved analysis showed that H I, Si III, and C IV roughly trace cold gas inside the Milky Way and Andromeda haloes, while O VI, O VII, and O VIII preferentially trace hot halo and intragroup gas, and that pressures outside the Milky Way halo are systematically higher toward the barycentre direction by a factor $\delta_{\rm RR}$3–$\delta_{\rm RR}$4 for $\delta_{\rm RR}$5 (Biaus et al., 19 Jan 2026).

The constrained Local Group also frames the stellar-halo studies. In six HESTIA galaxies, each system experiences one to four significant mergers with stellar mass ratio $\delta_{\rm RR}$6, most of them $\delta_{\rm RR}$7–$\delta_{\rm RR}$8 ago, and these events produce jumps in eccentricity and drops in $\delta_{\rm RR}$9 that reproduce Splash- and Plume-like kinematic features (Khoperskov et al., 2022). The in-situ component contributes about $$\delta_{\rm CR}(\mathbf{x})=\delta_{\rm RR}(\mathbf{x})+\sum_{i,j}\langle\delta(\mathbf{x})\,c_i\rangle\,[\langle c_i\,c_j\rangle]^{-1}(c_j-\tilde c_j),$$0–$$\delta_{\rm CR}(\mathbf{x})=\delta_{\rm RR}(\mathbf{x})+\sum_{i,j}\langle\delta(\mathbf{x})\,c_i\rangle\,[\langle c_i\,c_j\rangle]^{-1}(c_j-\tilde c_j),$$1 of the stellar halo mass within $$\delta_{\rm CR}(\mathbf{x})=\delta_{\rm RR}(\mathbf{x})+\sum_{i,j}\langle\delta(\mathbf{x})\,c_i\rangle\,[\langle c_i\,c_j\rangle]^{-1}(c_j-\tilde c_j),$$2 and $$\delta_{\rm CR}(\mathbf{x})=\delta_{\rm RR}(\mathbf{x})+\sum_{i,j}\langle\delta(\mathbf{x})\,c_i\rangle\,[\langle c_i\,c_j\rangle]^{-1}(c_j-\tilde c_j),$$3 beyond $$\delta_{\rm CR}(\mathbf{x})=\delta_{\rm RR}(\mathbf{x})+\sum_{i,j}\langle\delta(\mathbf{x})\,c_i\rangle\,[\langle c_i\,c_j\rangle]^{-1}(c_j-\tilde c_j),$$4 (Khoperskov et al., 2022). The accreted component occupies broad and overlapping structures in $$\delta_{\rm CR}(\mathbf{x})=\delta_{\rm RR}(\mathbf{x})+\sum_{i,j}\langle\delta(\mathbf{x})\,c_i\rangle\,[\langle c_i\,c_j\rangle]^{-1}(c_j-\tilde c_j),$$5–$$\delta_{\rm CR}(\mathbf{x})=\delta_{\rm RR}(\mathbf{x})+\sum_{i,j}\langle\delta(\mathbf{x})\,c_i\rangle\,[\langle c_i\,c_j\rangle]^{-1}(c_j-\tilde c_j),$$6, the UV plane, and action space, while a kinematic space based on $$\delta_{\rm CR}(\mathbf{x})=\delta_{\rm RR}(\mathbf{x})+\sum_{i,j}\langle\delta(\mathbf{x})\,c_i\rangle\,[\langle c_i\,c_j\rangle]^{-1}(c_j-\tilde c_j),$$7 and eccentricity separates debris more clearly (Khoperskov et al., 2022). The chemical analysis further shows that accreted debris are chemically distinct from surviving dwarf galaxies and that multi-element abundance patterns, together with stellar ages, are required to disentangle overlapping merger debris (Khoperskov et al., 2022).

Several other investigations rely directly on the constrained cosmography. Satellite accretion onto the Milky Way and Andromeda analogues is strongly aligned with the slowest-collapse axis $$\delta_{\rm CR}(\mathbf{x})=\delta_{\rm RR}(\mathbf{x})+\sum_{i,j}\langle\delta(\mathbf{x})\,c_i\rangle\,[\langle c_i\,c_j\rangle]^{-1}(c_j-\tilde c_j),$$8 of the tidal and shear tensors, especially for the early-infall population at $$\delta_{\rm CR}(\mathbf{x})=\delta_{\rm RR}(\mathbf{x})+\sum_{i,j}\langle\delta(\mathbf{x})\,c_i\rangle\,[\langle c_i\,c_j\rangle]^{-1}(c_j-\tilde c_j),$$9, and satellites can travel up to $\tilde c_j$0 relative to their parent halo before crossing $\tilde c_j$1 (Dupuy et al., 2022). The gravitational potential of Milky Way analogues is measurably non-stationary: at distances $\tilde c_j$2, ignoring the mass distribution outside the virial radius produces $\tilde c_j$3 errors in the potential quadrupole, and the spherical-harmonic coefficients vary significantly during the last $\tilde c_j$4 (Arakelyan et al., 2024). Radiative-transfer post-processing of the constrained Local Group shows that reionization proceeds in an inside-out manner, with the progenitors of the Milky Way and M31 reaching $\tilde c_j$5 ionization at $\tilde c_j$6–$\tilde c_j$7, earlier than the global midpoint at $\tilde c_j$8–$\tilde c_j$9, while external fronts contribute less than $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$00 of the ionized volume in the zoom region even in the least-suppressed model (Attard et al., 12 Sep 2025). In a Magellanic-analog dwarf pair, the more massive dwarf hosts a warm coronal envelope with $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$01, tidal interactions produce an H I stream of $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$02 and $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$03, and warm coronal gas is found to be ubiquitous in all halos with $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$04 in the HESTIA suite (Chisholm et al., 21 Apr 2025).

The same constrained simulations have also been used to study the inner Milky Way dark-matter morphology. In that application, the inner dark-matter halo is triaxial and boxy, with typical minor-to-major axis ratios at the $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$05 contour of $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$06–$$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$07 and $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$08–$$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$09, implying a non-spherical annihilation morphology for the Galactic Center excess (Muru et al., 8 Aug 2025). A plausible implication is that the constrained merger history and Local Group environment provide a dynamically specific context for predictions that would be obscured in unconstrained Milky Way analogues.

6. Interpretation, limitations, and role within Local-Volume cosmology

HESTIA is frequently used to argue that environmental effects are non-negligible even for apparently galaxy-scale observables. The suite was explicitly designed to preserve the large-scale tidal field, local filament/void environment, and Virgo-like surroundings while allowing different small-scale random modes between realizations (Arakelyan et al., 2024). This enables controlled comparisons between shared cosmography and realization-to-realization diversity.

At the same time, several analyses identify limitations. The 2022 Local Group kinematics study noted that the high-resolution HESTIA set lacks a Magellanic-Stream analogue at $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$10, that disentangling circumgalactic from intragroup absorbers remains challenging because they overlap in velocity and angular position, and that the study did not include detailed ionization modeling of the gas (Biaus et al., 2022). The stellar-halo studies emphasized that a six-galaxy sample limits statistical conclusions, that the resolution may miss low-mass substructure and fine-grained phase-space structure, and that broad abundance spreads in the adopted model complicate detailed chemical tagging (Khoperskov et al., 2022). The gas-column study likewise treated Milky Way contamination of M31 sightlines as a serious interpretive issue rather than a numerical inconsistency (Damle et al., 2022).

Within constrained-simulation methodology more broadly, the HESTIA program can be understood as the hydrodynamical extension of the RZA+CR strategy validated by Doumler et al. The essential advance is the backward Lagrangian shift of velocity constraints, which removes the $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$11 cosmic shift and enables recovery of $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$12-scale structures with $$\mathbf{v}(\mathbf{x}) = -H_0\,f\,\nabla \Phi(\mathbf{x}),$$13 positional accuracy in mock tests (Doumler et al., 2012). In HESTIA, this principle is combined with high-resolution zoom simulations and full baryonic physics so that the Local Group is modeled not as an isolated halo pair, but as a pair embedded in the observed cosmographic landscape (Libeskind et al., 2020).