Mira–Titan Universe Simulation Suite

• Mira–Titan is a cosmological simulation framework using gravity-only N-body runs and emulators to predict large-scale structure observables over an 8-dimensional parameter space.
• It employs a multi-fidelity design that combines high-resolution simulations, lower-resolution ensembles, and perturbation theory to achieve 2–3% precision in matter power spectrum and halo mass function predictions.
• The suite serves both as a computational asset for dark energy surveys and as a testbed for advanced statistical methods, including Bayesian deep Gaussian process modeling.

Searching arXiv for recent and foundational papers on the Mira-Titan Universe simulation suite and related emulation work. The Mira–Titan Universe Simulation Suite is a cosmological simulation and emulation program built to provide precision predictions for nonlinear large-scale structure observables over an 8-dimensional cosmological parameter space that extends beyond standard $\Lambda$CDM to include massive neutrinos and dynamical dark energy. Across the literature, the suite is described as a carefully designed ensemble of gravity-only $N$-body simulations run with HACC on the Mira and Titan supercomputers, with a nested, space-filling design intended to support statistically controlled emulators for the matter power spectrum, halo mass function, and later galaxy clustering and lensing observables (Heitmann et al., 2015). Subsequent work presents both the original multi-fidelity power-spectrum construction and its final high-precision realization, as well as newer statistical front ends that treat the suite’s outputs as correlated functional data and explicitly quantify uncertainty at the single-cosmology level (Lawrence et al., 2017, Moran et al., 2022, Walsh et al., 24 Jul 2025).

1. Historical development and scientific purpose

The suite was introduced as a next-generation precision cosmology program for dark energy surveys, with explicit targets including weak gravitational lensing, baryon acoustic oscillations, galaxy cluster abundance, redshift-space distortions, and galaxy clustering (Heitmann et al., 2015). Its central methodological premise is that full high-resolution simulations cannot be run densely across cosmological parameter space, whereas a limited set of carefully placed simulations can be used to train Gaussian-process-based emulators that deliver fast predictions with controlled and incrementally improvable error (Heitmann et al., 2015).

Within the Mira–Titan series, the program matured through a sequence of emulator products. The initial design paper established the 8-dimensional parameter space and nested lattice strategy (Heitmann et al., 2015). The first matter power spectrum emulator used the first 36 cosmologies and combined perturbation theory, medium-resolution simulations, and high-resolution $N$-body runs, achieving approximately $4\%$ accuracy over $0 \le z \le 2$ and $k \le 5\,\mathrm{Mpc}^{-1}$ (Lawrence et al., 2017). The later “Mira-Titan Universe IV” paper presented the final high-precision matter power spectrum emulator using 111 cosmological simulations, 1776 lower-resolution simulations, and TimeRG perturbation theory results, with predictions at the two to three percent level over the full design domain (Moran et al., 2022). The halo mass function emulator extended the same simulation backbone to group and cluster scales, obtaining typical precision better than $2\%$ for $10^{13}$–$10^{14}\,M_\odot/h$ at $z \lesssim 1$ and $N$0 at $N$1 (Bocquet et al., 2020).

A later statistical reinterpretation treats the suite not merely as a source of emulator training targets, but as a structured collection of correlated, multi-fidelity functional outputs. In that formulation, the Bayesian deep Gaussian process hierarchical model is explicitly built to exploit the strengths and limitations of Mira–Titan power spectrum outputs, first refining raw simulation curves into uncertainty-quantified latent spectra and then training a separate emulator across cosmologies (Walsh et al., 24 Jul 2025). This suggests that the suite has functioned both as a computational asset and as a testbed for increasingly sophisticated statistical methodology.

2. Cosmological parameter space and design structure

The suite spans an 8-dimensional parameter space comprising

$N$2

or equivalently, in later summaries, matter density, baryon density, amplitude of density fluctuations, dimensionless Hubble parameter, scalar spectral index, dark energy equation of state parameters, and neutrino density (Moran et al., 2022, Walsh et al., 24 Jul 2025). The final matter power spectrum paper gives the parameter ranges as

$N$3

together with the joint constraint

$N$4

while assuming $N$5 and fixing $N$6 (Moran et al., 2022).

Dark energy is modeled with the CPL form

$N$7

and the suite is explicitly constructed to cover both massless and massive neutrino cosmologies as well as dynamical dark energy (Heitmann et al., 2015, Moran et al., 2022). In the simulations, dark energy is treated as a smooth component affecting the background expansion, while neutrinos are handled approximately through linear theory and background evolution rather than particle-based nonlinear clustering (Moran et al., 2022, Bocquet et al., 2020).

The design is a defining feature. The original papers describe a space-filling lattice design with nested refinement, allowing staged improvement of emulator fidelity without discarding earlier simulations (Heitmann et al., 2015). In the final 111-model realization, the sequence is summarized as 26 initial lattice-based models, plus 10 additional massless-neutrino models chosen via a symmetric Latin hypercube, then 29 new models, and finally 46 more, reaching 111 total (Moran et al., 2022). The 2025 statistical treatment uses 117 cosmologies from the Mira–Titan dataset, split into 111 for training and 6 held out for testing and benchmarking at $N$8 (Walsh et al., 24 Jul 2025).

This staged design has two consequences visible across the literature. First, it supports progressively refined emulators, from the 36-model intermediate-accuracy realization to the final 111-model products (Lawrence et al., 2017, Moran et al., 2022). Second, it creates a structured training set in which later papers can compare alternative statistical front ends while holding fixed the underlying cosmological design (Walsh et al., 24 Jul 2025).

3. Simulation campaign and numerical realization

The suite consists of gravity-only HACC simulations executed on Mira and Titan, with the final matter power spectrum paper describing one high-resolution realization per cosmology, supplemented by multiple lower-resolution realizations and perturbation theory (Moran et al., 2022). For each of the 111 cosmological models in the final design, the numerical setup includes a high-resolution simulation with volume $N$9 and $N$0 CDM+baryon particles, and a set of lower-resolution simulations with volume $N$1; across the full suite there are 1776 lower-resolution runs, and 16 lower-resolution realizations are used in the smoothing stage for each cosmology/redshift combination (Moran et al., 2022).

Earlier power-spectrum work used essentially the same multi-fidelity structure in an intermediate design with 36 cosmologies: one high-resolution run with $N$2 particles in a $N$3 box, sixteen medium-resolution particle–mesh realizations with $N$4 particles on a $N$5 mesh in a $N$6 box, and TimeRG perturbation theory on large scales (Lawrence et al., 2017). The original design study had already identified this matching strategy—perturbation theory at low $N$7, multiple lower-resolution realizations at intermediate $N$8, and one high-resolution run at high $N$9—as the appropriate route to percent-level calibration over $4\%$0 (Heitmann et al., 2015).

The final power-spectrum emulator covers eight redshifts,

$4\%$1

and is valid over $4\%$2 and up to $4\%$3 (Moran et al., 2022). By contrast, the 2025 Bayesian deep Gaussian process case study fixes redshift to $4\%$4 and focuses on the internal structure of the multi-fidelity power spectrum outputs for each cosmology (Walsh et al., 24 Jul 2025).

The suite is embedded in a broader HACC ecosystem that includes very large single-cosmology simulations such as Outer Rim and Last Journey. Those simulations are not themselves members of the parameter-space ensemble, but they use the same code base, analysis tools, and emulator infrastructure, and are explicitly cross-validated against Mira–Titan emulators (Heitmann et al., 2019, Heitmann et al., 2020). A plausible implication is that Mira–Titan should be understood as the ensemble arm of a larger HACC-based program in which survey-scale flagship runs provide complementary validation and mock-generation capability.

4. Multi-fidelity matter power spectra and emulation workflow

The matter power spectrum is the suite’s most developed observable. In the final high-precision realization, the emulator is built by combining TimeRG perturbation theory, 16 lower-resolution simulations, and one high-resolution simulation per cosmology, aligned across overlapping $4\%$5-ranges and transformed into the emulation space

$4\%$6

before smoothing and emulation (Moran et al., 2022). The final data combination uses perturbation theory for $4\%$7, both lower-resolution and high-resolution data for $4\%$8, and high-resolution data only for $4\%$9 (Moran et al., 2022).

The original 36-model emulator smoothed the combined spectra with a two-layer process convolution model and then represented the resulting $0 \le z \le 2$0 at each redshift through principal components whose weights were modeled as Gaussian processes over the 8-dimensional cosmological parameter space (Lawrence et al., 2017). The final 111-model emulator generalized this approach by stacking the 8 redshift spectra per cosmology, performing PCA on vectors of length $0 \le z \le 2$1, retaining $0 \le z \le 2$2 principal components, and emulating the corresponding weights with the SEPIA package (Moran et al., 2022).

The 2025 reinterpretation uses the same multi-fidelity scientific structure but introduces a different front end. For each cosmology at $0 \le z \le 2$3, the suite yields 18 power spectrum curves: one perturbation theory spectrum $0 \le z \le 2$4, sixteen low-resolution $0 \le z \le 2$5-body spectra $0 \le z \le 2$6, and one high-resolution $0 \le z \le 2$7-body spectrum $0 \le z \le 2$8, all evaluated on a common grid of $0 \le z \le 2$9 wavenumbers in $k \le 5\,\mathrm{Mpc}^{-1}$0 (Walsh et al., 24 Jul 2025). The paper defines a precision-weighted average

$k \le 5\,\mathrm{Mpc}^{-1}$1

with $k \le 5\,\mathrm{Mpc}^{-1}$2 diagonal precision matrices and $k \le 5\,\mathrm{Mpc}^{-1}$3 the mean of the 16 low-resolution runs (Walsh et al., 24 Jul 2025). The crucial observation is that the low-resolution runs are not independent pointwise noise but smooth correlated draws around a latent true spectrum.

On that basis, the paper models the underlying true spectrum $k \le 5\,\mathrm{Mpc}^{-1}$4 through a three-layer deep Gaussian process with functional correlated outputs: $k \le 5\,\mathrm{Mpc}^{-1}$5 where $k \le 5\,\mathrm{Mpc}^{-1}$6, $k \le 5\,\mathrm{Mpc}^{-1}$7 is a latent monotone warp, and the covariance $k \le 5\,\mathrm{Mpc}^{-1}$8 encodes both wavenumber-dependent fidelity and dense low-resolution correlation structure (Walsh et al., 24 Jul 2025). The warp is used to capture nonstationarity, especially around the BAO region $k \le 5\,\mathrm{Mpc}^{-1}$9 (Walsh et al., 24 Jul 2025).

After posterior inference for each training cosmology, the posterior mean spectra $2\%$0 are compressed with a PCA basis using $2\%$1 principal components, and the PC weights are emulated across the 8-dimensional cosmological parameter vector with scalar Gaussian processes using a power exponential correlation

$2\%$2

with $2\%$3 fixed and parameters estimated via maximum likelihood with the GPfit R package (Walsh et al., 24 Jul 2025). This produces a back-end emulator structurally similar to earlier Mira–Titan emulators, but trained on uncertainty-refined latent spectra rather than directly smoothed simulation outputs.

5. Halo mass function, galaxy observables, and extension beyond $2\%$4

The suite was designed from the outset to support emulators beyond the matter power spectrum, including the halo mass function and eventually halo concentration–mass relations, galaxy clustering, and redshift-space distortions (Heitmann et al., 2015). The halo mass function emulator uses the completed 111-model suite, extracting $2\%$5 halo catalogs in the redshift range $2\%$6 and for masses $2\%$7 (Bocquet et al., 2020). Halos are identified with a friends-of-friends stage using linking length $2\%$8, followed by spherical overdensity masses grown around the FOF potential minimum (Bocquet et al., 2020).

Instead of emulating a universal form $2\%$9, that work fits smooth piecewise second-order polynomials to $10^{13}$0, imposes continuity and derivative constraints between segments, compresses the resulting mass functions with PCA, and emulates the leading four principal component weights with Gaussian processes (Bocquet et al., 2020). The paper explicitly shows that traditional universal fitting forms can be biased at up to $10^{13}$1 at $10^{13}$2 for $10^{13}$3, whereas the Mira–Titan emulator maintains much smaller errors across the design space (Bocquet et al., 2020).

A further extension uses the 111-simulation suite as the basis for emulators of the projected and redshift-space galaxy correlation functions and excess surface density from galaxy-galaxy lensing, all within HOD modeling and all varying over the same eight cosmological parameters, including neutrino mass and dynamical dark energy (Kwan et al., 2023). In that work, each simulation has comoving box length $10^{13}$4, $10^{13}$5 particles, force resolution $10^{13}$6, and halos identified with FOF linking length $10^{13}$7 and a minimum of 40 particles per halo (Kwan et al., 2023). The emulators are built at a single snapshot $10^{13}$8, chosen to approximate the BOSS DR12 CMASS redshift distribution (Kwan et al., 2023).

The galaxy-clustering emulators combine an extended Zheng et al. HOD with velocity bias, redshift-space mapping, and a Kronecker-structured Gaussian process that exploits separability between cosmology and HOD design (Kwan et al., 2023). The resulting models are reported to be sufficiently accurate for analysis of the BOSS DR12 CMASS sample over $10^{13}$9, and mock tests yield $10^{14}\,M_\odot/h$0 constraints on the growth rate and $10^{14}\,M_\odot/h$1 on $10^{14}\,M_\odot/h$2 for a CMASS-like sample using only the measurements covered by the emulator, improving to a $10^{14}\,M_\odot/h$3 growth-rate measurement with an external $10^{14}\,M_\odot/h$4 prior (Kwan et al., 2023).

These developments show that Mira–Titan is not simply a power-spectrum suite. It is a common simulation substrate for multiple nonlinear observables, each modeled with a different statistical representation but sharing the same parameter-space coverage, halo definitions, and HACC numerical backbone.

6. Accuracy, comparisons, limitations, and broader significance

The final high-precision matter power spectrum emulator provides predictions at the two to three percent level over the full 8-dimensional domain, with out-of-sample test cosmologies reproduced to better than $10^{14}\,M_\odot/h$5 across $10^{14}\,M_\odot/h$6 and $10^{14}\,M_\odot/h$7 in all cases except one, which reaches $10^{14}\,M_\odot/h$8; for the WMAP7 $10^{14}\,M_\odot/h$9CDM test cosmology M000, accuracy is $z \lesssim 1$0 at all $z \lesssim 1$1 and redshifts (Moran et al., 2022). That paper also reports significant improvement over the earlier 36-model emulator, whose errors were $z \lesssim 1$2–$z \lesssim 1$3, and shows the final emulator outperforming Halofit, EuclidEmulator2, and HMCODE-2020 in the reported comparisons (Moran et al., 2022).

For the halo mass function emulator, the typical precision is better than $z \lesssim 1$4 for $z \lesssim 1$5–$z \lesssim 1$6 at $z \lesssim 1$7, with degradation toward the high-mass tail where halo counts are intrinsically sparse (Bocquet et al., 2020). For galaxy clustering and galaxy-galaxy lensing, the emulators reach roughly $z \lesssim 1$8 accuracy for $z \lesssim 1$9 and $N$00, and $N$01–$N$02 for redshift-space multipoles over the ranges analyzed (Kwan et al., 2023).

The 2025 Bayesian deep Gaussian process treatment does not replace the final public matter power spectrum emulator; rather, it offers an alternative statistical front end tailored to the correlated functional structure of the raw Mira–Titan outputs (Walsh et al., 24 Jul 2025). Its comparison is therefore not primarily against simulation truth for Mira–Titan, which is unavailable, but against Cosmic Emu and against in-sample DGP.FCO fits to held-out cosmologies (Walsh et al., 24 Jul 2025). On six held-out Mira–Titan cosmologies, DGP.FCO+PC achieves the lowest mean squared error on scales where perturbation theory and low-resolution runs are both approximately unbiased, while Cosmic Emu has a slight advantage where only the high-resolution run is unbiased; overall, DGP.FCO+PC has lower MSE at $N$03 of the $N$04-values considered (Walsh et al., 24 Jul 2025). This suggests methodological competitiveness rather than definitive replacement.

Several limitations recur across the literature. The suite is gravity-only, so baryonic effects are absent and must be modeled separately, especially for $N$05 or for low-mass halo statistics (Moran et al., 2022, Bocquet et al., 2020). Massive neutrinos are treated approximately, with linear neutrino clustering and no neutrino particles in the $N$06-body evolution, although the cited papers argue that this is adequate for the target mass and scale ranges (Moran et al., 2022, Bocquet et al., 2020). Predictions are valid only within the design hypercube, and extrapolation outside it is explicitly discouraged (Lawrence et al., 2017, Bocquet et al., 2020). Some later extensions also operate at only a single redshift, such as $N$07 for the galaxy-clustering emulators and $N$08 for the DGP.FCO case study (Kwan et al., 2023, Walsh et al., 24 Jul 2025).

At the same time, the suite’s broader significance lies in its combination of numerical design and emulator-oriented statistical structure. It provides a shared 8-dimensional cosmological response surface on which multiple observables can be modeled consistently, from $N$09 to the halo mass function to HOD-based clustering and lensing (Heitmann et al., 2015, Moran et al., 2022, Bocquet et al., 2020, Kwan et al., 2023). Later work further suggests that when a simulation suite delivers multiple correlated functional runs per input, with well-characterized fidelity structure, one can build statistically richer front ends than simple smoothing or diagonal-noise Gaussian processes (Walsh et al., 24 Jul 2025). In that sense, Mira–Titan is both a simulation campaign and a reference architecture for emulator construction in nonlinear cosmology.