Initial Conditions For N-body Simulations Research Articles

Bayesian inference of initial conditions from non-linear cosmic structures using field-level emulators

ABSTRACT Analysing next-generation cosmological data requires balancing accurate modelling of non-linear gravitational structure formation and computational demands. We propose a solution by introducing a machine learning-based field-level emulator, within the Hamiltonian Monte Carlo-based Bayesian Origin Reconstruction from Galaxies (BORG) inference algorithm. Built on a V-net neural network architecture, the emulator enhances the predictions by first-order Lagrangian perturbation theory to be accurately aligned with full N-body simulations while significantly reducing evaluation time. We test its incorporation in BORG for sampling cosmic initial conditions using mock data based on non-linear large-scale structures from N-body simulations and Gaussian noise. The method efficiently and accurately explores the high-dimensional parameter space of initial conditions, fully extracting the cross-correlation information of the data field binned at a resolution of $1.95\,h^{-1}$ Mpc. Percent-level agreement with the ground truth in the power spectrum and bispectrum is achieved up to the Nyquist frequency $k_\mathrm{N} \approx 2.79h \,\, \mathrm{Mpc}^{-1}$. Posterior resimulations – using the inferred initial conditions for N-body simulations – show that the recovery of information in the initial conditions is sufficient to accurately reproduce halo properties. In particular, we show highly accurate $M_{200\mathrm{c}}$ halo mass function and stacked density profiles of haloes in different mass bins $[0.853,16]\times 10^{14}\,{\rm M}_{\odot }\,h^{-1}$. As all available cross-correlation information is extracted, we acknowledge that limitations in recovering the initial conditions stem from the noise level and data grid resolution. This is promising as it underscores the significance of accurate non-linear modelling, indicating the potential for extracting additional information at smaller scales.

Algorithm to Produce a Density Field with Given Two, Three, and Four-Point Correlation Functions

Abstract Here we show how to produce a 3D density field with a given set of higher-order correlation functions. Our algorithm enables producing any desired two-point, three-point, and four-point functions, including odd-parity for the latter. We note that this algorithm produces the desired correlations about a set of “primary” points, matched to how the spherical-harmonic-based algorithms ENCORE and CADENZA measure them. These “primary points” must be used as those around which the correlation functions are measured. We also generalize the algorithm to i) N-point correlations with N > 4, ii) dimensions other than 3, and iii) beyond scalar quantities. This algorithm should find use in verifying analysis pipelines for higher-order statistics in upcoming galaxy redshift surveys such as DESI, Euclid, Roman, and Spherex, as well as intensity mapping. In particular it may be helpful in searches for parity violation in the 4PCF of these samples, for which producing initial conditions for N-body simulations is both costly and highly model-dependent at present, and so alternative methods such as that developed here are desirable.

Comparisons of the core and mantle compositions of earth analogs from different terrestrial planet formation scenarios

The chemical compositions of Earth's core and mantle provide insight into the processes that led to their formation. N-body simulations, on the other hand, generally do not contain chemical information, and seek to only reproduce the masses and orbits of the terrestrial planets. These simulations can be grouped into four potentially viable scenarios of Solar System formation (Classical, Annulus, Grand Tack, and Early Instability) for which we compile a total of 433 N-body simulations. We relate the outputs of these simulations to the chemistry of Earth's core and mantle using a melt-scaling law combined with a multi-stage model of core formation. We find the compositions of Earth analogs to be largely governed by the fraction of equilibrating embryo cores (kcore_emb) and the initial embryo masses in N-body simulations, rather than the simulation type, where higher values of kcore_emb and larger initial embryo masses correspond to higher concentrations of Ni, Co, Mo, and W in Earth analog mantles and higher concentrations of Si and O in Earth analog cores. As a result, larger initial embryo masses require smaller values of kcore_emb to match Earth's mantle composition. On the other hand, compositions of Earth analog cores are sensitive to the temperatures of equilibration and fO2 of accreting material. Simulation type may be important when considering magma ocean lifetimes, where Grand Tack simulations have the largest amounts of material accreted after the last giant impact. However, we cannot rule out any accretion scenarios or initial embryo masses due to the sensitivity of Earth's mantle composition to different parameters and the stochastic nature of N-body simulations. We use our compiled simulations to explore the relationship between initial embryo masses and the melting history of Earth analogs, where the complex interplay between the timing between impacts, magma ocean lifetimes, and volatile delivery could affect the compositions of Earth analogs formed from different simulation types. Comparing the last embryo impacts experienced by Earth analogs to specific Moon-forming scenarios, we find the characteristics of the Moon-forming impact are dependent on the initial conditions in N-body simulations where larger initial embryo masses promote larger and slower Moon-forming impactors. Mars-sized initial embryos are most consistent with the canonical hit-and-run scenario onto a solid mantle. Our results suggest that constraining the fraction of equilibrating impactor core (kcore) and the initial embryo masses in N-body simulations could be significant for understanding both Earth's accretion history and characteristics of the Moon-forming impact.

Simple lessons from complex learning: what a neural network model learns about cosmic structure formation.

We train a neural network model to predict the full phase space evolution of cosmological N-body simulations. Its success implies that the neural network model is accurately approximating the Green's function expansion that relates the initial conditions of the simulations to its outcome at later times in the deeply nonlinear regime. We test the accuracy of this approximation by assessing its performance on well-understood simple cases that have either known exact solutions or well-understood expansions. These scenarios include spherical configurations, isolated plane waves, and two interacting plane waves: initial conditions that are very different from the Gaussian random fields used for training. We find our model generalizes well to these well-understood scenarios, demonstrating that the networks have inferred general physical principles and learned the nonlinear mode couplings from the complex, random Gaussian training data. These tests also provide a useful diagnostic for finding the model's strengths and weaknesses, and identifying strategies for model improvement. We also test the model on initial conditions that contain only transverse modes, a family of modes that differ not only in their phases but also in their evolution from the longitudinal growing modes used in the training set. When the network encounters these initial conditions that are orthogonal to the training set, the model fails completely. In addition to these simple configurations, we evaluate the model's predictions for the density, displacement, and momentum power spectra with standard initial conditions for N-body simulations. We compare these summary statistics against N-body results and an approximate, fast simulation method called COLA (COmoving Lagrangian Acceleration). Our model achieves percent level accuracy at nonlinear scales of , representing a significant improvement over COLA.

From hydrodynamics to N-body simulations of star clusters: mergers and rotation

ABSTRACT We present a new method to obtain more realistic initial conditions for N-body simulations of young star clusters. We start from the outputs of hydrodynamical simulations of molecular cloud collapse, in which star formation is modelled with sink particles. In our approach, we instantaneously remove gas from these hydrodynamical simulation outputs to mock the end of the gas-embedded phase, induced by stellar feedback. We then enforce a realistic initial mass function by splitting or joining the sink particles based on their mass and position. Such initial conditions contain more consistent information on the spatial distribution and the kinematical and dynamical states of young star clusters, which are fundamental to properly study these systems. For example, by applying our method to a set of previously run hydrodynamical simulations, we found that the early evolution of young star clusters is affected by gas removal and by the early dry merging of sub-structures. This early evolution can either quickly erase the rotation acquired by our (sub-)clusters in their embedded phase or ‘fuel’ it by feeding of angular momentum by sub-structure mergers, before two-body relaxation acts on longer time-scales

Predicting dark matter halo formation in N-body simulations with deep regression networks

ABSTRACT Dark matter haloes play a fundamental role in cosmological structure formation. The most common approach to model their assembly mechanisms is through N-body simulations. In this work, we present an innovative pathway to predict dark matter halo formation from the initial density field using a Deep Learning algorithm. We implement and train a Deep Convolutional Neural Network to solve the task of retrieving Lagrangian patches from which dark matter haloes will condense. The volumetric multilabel classification task is turned into a regression problem by means of the Euclidean distance transformation. The network is complemented by an adaptive version of the watershed algorithm to form the entire protohalo identification pipeline. We show that splitting the segmentation problem into two distinct subtasks allows for training smaller and faster networks, while the predictive power of the pipeline remains the same. The model is trained on synthetic data derived from a single full N-body simulation and achieves deviations of ∼10 per cent when reconstructing the dark matter halo mass function at z = 0. This approach represents a promising framework for learning highly non-linear relations in the primordial density field. As a practical application, our method can be used to produce mock dark matter halo catalogues directly from the initial conditions of N-body simulations.

AGAMA: action-based galaxy modelling architecture

Agama is a publicly available software library for a broad range of applications in the field of stellar dynamics. It provides methods for computing the gravitational potential of arbitrary analytic density profiles or N-body models; orbit integration and analysis; transformations between position/velocity and action/angle variables; distribution functions expressed in terms of actions and their moments; iterative construction of self-consistent multicomponent galaxy models. Applications include the inference about the structure of Milky Way or other galaxies from observations of stellar kinematics; preparation of equilibrium initial conditions for N-body simulations; analysis of snapshots from simulations. The library is written in C++, provides a Python interface, and can be coupled to other stellar-dynamical software: Amuse, Galpy and Nemo.

Formation of solar system analogues – I. Looking for initial conditions through a population synthesis analysis

Population synthesis models of planetary systems developed during the last $\sim$15 years could reproduce several of the observables of the exoplanet population, and also allowed to constrain planetary formation models. We present our planet formation model, which calculates the evolution of a planetary system during the gaseous phase. The code incorporates relevant physical phenomena for the formation of a planetary system, like photoevaporation, planet migration, gas accretion, water delivery in embryos and planetesimals, a detailed study of the orbital evolution of the planetesimal population, and the treatment of the fusion between embryos, considering their atmospheres. The main goal of this work, unlike other works of planetary population synthesis, is to find suitable scenarios and physical parameters of the disc to form solar system analogs. We are specially interested in the final planet distributions, and in the final surface density, eccentricity and inclination profiles for the planetesimal population. These final distributions will be used as initial conditions for N-body simulations, to study the post-oligarchic formation in a second work. We then consider different formation scenarios, with different planetesimal sizes and different type I migration rates. We find that solar system analogs are favored in massive discs, with low type I migration rates, and small planetesimal sizes. Besides, those rocky planets within their habitables zones are dry when discs dissipate. At last, the final configurations of solar system analogs include information about the mass and semimajor-axis of the planets, water contents, and the properties of the planetesimal remnants.

SPH/N-Body simulations of small ([formula omitted]) asteroidal breakups and improved parametric relations for Monte–Carlo collisional models

We report on our study of asteroidal breakups, i.e. fragmentations of targets, subsequent gravitational reaccumulation and formation of small asteroid families. We focused on parent bodies with diameters Dpb=10km. Simulations were performed with a smoothed-particle hydrodynamics (SPH) code combined with an efficient N-body integrator. We assumed various projectile sizes, impact velocities and impact angles (125 runs in total). Resulting size-frequency distributions are significantly different from scaled-down simulations with Dpb=100km targets (Durda et al., 2007). We derive new parametric relations describing fragment distributions, suitable for Monte-Carlo collisional models. We also characterize velocity fields and angular distributions of fragments, which can be used as initial conditions for N-body simulations of small asteroid families. Finally, we discuss a number of uncertainties related to SPH simulations.

Relativistic initial conditions for N-body simulations

Initial conditions for (Newtonian) cosmological N-body simulations are usually set by re-scaling the present-day power spectrum obtained from linear (relativistic) Boltzmann codes to the desired initial redshift of the simulation. This back-scaling method can account for the effect of inhomogeneous residual thermal radiation at early times, which is absent in the Newtonian simulations. We analyse this procedure from a fully relativistic perspective, employing the recently-proposed Newtonian motion gauge framework. We find that N-body simulations for ΛCDM cosmology starting from back-scaled initial conditions can be self-consistently embedded in a relativistic space-time with first-order metric potentials calculated using a linear Boltzmann code. This space-time coincides with a simple ``N-body gauge'' for z < 50 for all observable modes. Care must be taken, however, when simulating non-standard cosmologies. As an example, we analyse the back-scaling method in a cosmology with decaying dark matter, and show that metric perturbations become large at early times in the back-scaling approach, indicating a breakdown of the perturbative description. We suggest a suitable ``forwards approach" for such cases.

Squeezing the halo bispectrum: a test of bias models

We study the halo-matter cross bispectrum in the presence of primordial non-Gaussianity of the local type. We restrict ourselves to the squeezed limit, for which the calculation are straightforward, and perform the measurements in the initial conditions of N-body simulations, to mitigate the contamination induced by nonlinear gravitational evolution. Interestingly, the halo-matter cross bispectrum is not trivial even in this simple limit as it is strongly sensitive to the scale-dependence of the quadratic and third-order halo bias. Therefore, it can be used to test biasing prescriptions. We consider three different prescription for halo clustering: excursion set peaks (ESP), local bias and a model in which the halo bias parameters are explicitly derived from a peak-background split. In all cases, the model parameters are fully constrained with statistics other than the cross bispectrum. We measure the cross bispectrum involving one halo fluctuation field and two mass overdensity fields for various halo masses and collapse redshifts. We find that the ESP is in reasonably good agreement with the numerical data, while the other alternatives we consider fail in various cases. This suggests that the scale-dependence of halo bias also is a crucial ingredient to the squeezed limit of the halo bispectrum.

The Impact of Initial Conditions in N-Body Simulations of Debris Discs

AbstractNumerical simulations are a crucial tool to understand the relationship between debris discs and planetary companions. As debris disc observations are now reaching unprecedented levels of precision over a wide range of wavelengths, an appropriate level of accuracy and consistency is required in numerical simulations to confidently interpret this new generation of observations. However, simulations throughout the literature have been conducted with various initial conditions often with little or no justification. In this paper, we aim to study the dependence on the initial conditions of N-body simulations modelling the interaction between a massive and eccentric planet on an exterior debris disc. To achieve this, we first classify three broad approaches used in the literature and provide some physical context for when each category should be used. We then run a series of N-body simulations, that include radiation forces acting on small grains, with varying initial conditions across the three categories. We test the influence of the initial parent body belt width, eccentricity, and alignment with the planet on the resulting debris disc structure and compare the final peak emission location, disc width and offset of synthetic disc images produced with a radiative transfer code. We also track the evolution of the forced eccentricity of the dust grains induced by the planet, as well as resonance dust trapping. We find that an initially broad parent body belt always results in a broader debris disc than an initially narrow parent body belt. While simulations with a parent body belt with low initial eccentricity (e ~ 0) and high initial eccentricity (0 &lt; e &lt; 0.3) resulted in similar broad discs, we find that purely secular forced initial conditions, where the initial disc eccentricity is set to the forced value and the disc is aligned with the planet, always result in a narrower disc. We conclude that broad debris discs can be modelled by using either a dynamically cold or dynamically warm parent belt, while in contrast eccentric narrow debris rings are reproduced using a secularly forced parent body belt.

NON-GAUSSIAN INITIAL CONDITIONS IN ΛCDM: NEWTONIAN, RELATIVISTIC, AND PRIMORDIAL CONTRIBUTIONS

The goal of the present paper is to set initial conditions for structure formation at non-linear order, consistent with general relativity, while also allowing for primordial non-Gaussianity. We use the non-linear continuity and Raychauduri equations, which together with the non-linear energy constraint determine the evolution of the matter density fluctuation in general relativity. We solve this equations at first and second order in a perturbative expansion, recovering and extending previous results derived in the matter-dominated limit and in the Newtonian regime. We present a second-order solution for the comoving density contrast in a $\Lambda$CDM universe, identifying non-linear contributions coming from the Newtonian growing mode, primordial non-Gaussianity and intrinsic non-Gaussianity, due to the essential non-linearity of the relativistic constraint equations. We discuss the application of these results to initial conditions in N-body simulations, showing that relativistic corrections mimic a non-zero non-linear parameter $f_{\rm NL}$.

Cosmological structure formation with augmented Lagrangian perturbation theory

AbstractWe present a new fast and efficient approach to model structure formation with augmented Lagrangian perturbation theory (ALPT). Our method is based on splitting the displacement field into a long- and a short-range component. The long-range component is computed by second-order LPT (2LPT). This approximation contains a tidal non-local and non-linear term. Unfortunately, 2LPT fails on small scales due to severe shell crossing and a crude quadratic behaviour in the low-density regime. The spherical collapse (SC) approximation has been recently reported to correct for both effects by adding an ideal collapse truncation. However, this approach fails to reproduce the structures on large scales where it is significantly less correlated with the N-body result than 2LPT or linear LPT (the Zel'dovich approximation). We propose to combine both approximations using for the short-range displacement field the SC solution. A Gaussian filter with a smoothing radius rS is used to separate between both regimes. We use the result of 25 dark-matter-only N-body simulations to benchmark at z = 0 the different approximations: first-, second-, third-order LPT, SC and our novel combined ALPT model. This comparison demonstrates that our method improves previous approximations at all scales showing ∼25 and ∼75 per cent higher correlation than 2LPT with the N-body solution at k = 1 and 2 h Mpc−1, respectively. We conduct a parameter study to determine the optimal range of smoothing radii and find that the maximum correlation is achieved with rS = 4–5 h−1 Mpc. This structure formation approach could be used for various purposes, such as setting-up initial conditions for N-body simulations, generating mock galaxy catalogues, cosmic web analysis or for reconstructions of the primordial density fluctuations.

Three-integral multicomponent dynamical models and simulations of the nuclear star cluster in NGC 4244

Adaptive optics observations of the flattened nuclear star cluster in the nearby edge-on spiral galaxy NGC 4244 using the Gemini Near-Infrared Integral Field Spectrograph (NIFS) have revealed clear rotation. Using these kinematics plus 2MASS photometry we construct a series of axisymmetric two-component particle dynamical models with our improved version of NMAGIC, a flexible Chi^2-made-to-measure code. The models consist of a nuclear cluster disc embedded within a spheroidal particle population. We find a mass for the nuclear star cluster of M=1.6^+0.5_-0.2 x 10^7 M_sun within ~42.4 pc (2"). We also explore the presence of an intermediate mass black hole and show that models with a black hole as massive as M_bh = 5.0 x 10^5 M_sun are consistent with the available data. Regardless of whether a black hole is present or not, the nuclear cluster is vertically anisotropic (beta_z < 0), as was found with earlier two-integral models. We then use the models as initial conditions for N-body simulations. These simulations show that the nuclear star cluster is stable against non-axisymmetric perturbations. We also explore the effect of the nuclear cluster accreting star clusters at various inclinations. Accretion of a star cluster with mass 13% that of the nuclear cluster is already enough to destroy the vertical anisotropy, regardless of orbital inclination.

N-body simulations with generic non-Gaussian initial conditions I: power spectrum and halo mass function

We address the issue of setting up generic non-Gaussian initial conditions for N-body simulations. We consider inflationary-motivated primordial non-Gaussianity where the perturbations in the Bardeen potential are given by a dominant Gaussian part plus a non-Gaussian part specified by its bispectrum. The approach we explore here is suitable for any bispectrum, i.e. it does not have to be of the so-called separable or factorizable form. The procedure of generating a non-Gaussian field with a given bispectrum (and a given power spectrum for the Gaussian component) is not univocal, and care must be taken so that higher-order corrections do not leave a too large signature on the power spectrum. This is so far a limiting factor of our approach.We then run N-body simulations for the most popular inflationary-motivated non-Gaussian shapes. The halo mass function and the non-linear power spectrum agree with theoretical analytical approximations proposed in the literature, even if they were so far developed and tested only for a particular shape (the local one). We plan to make the simulations outputs available to the community via the non-Gaussian simulations comparison project web site http://icc.ub.edu/~liciaverde/NGSCP.html.

An iterative method for constructing equilibrium phase models of stellar systems

We present a new method for constructing equilibrium phase models for stellar systems, which we call the iterative method. It relies on constrained, or guided evolution, so that the equilibrium solution has a number of desired parameters and/or constraints. This method is very powerful, to a large extent due to its simplicity. It can be used for mass distributions with an arbitrary geometry and a large variety of kinematical constraints. We present several examples illustrating it. Applications of this method include the creation of initial conditions for N-body simulations and the modelling of galaxies from their photometric and kinematic observations.

Transients from initial conditions based on Lagrangian perturbation theory inN-body simulations

We explore the initial conditions for cosmologicalN-body simulations suitable for calculating the skewness and kurtosis of thedensity field. In general, the initial conditions based on the perturbation theoryprovide incorrect second-order and higher-order growth. These errors impliedby the use of the perturbation theory to set up the initial conditions inN-body simulations are called transients. Unless these transients are completely suppressed comparedwith the dominant growing mode, we cannot reproduce the correct evolution of cumulants withorders higher than two, even though there is no problem with the numerical scheme. Weinvestigate the impact of transients on the observable statistical quantities by performingN-body simulations with initial conditions based on Lagrangian perturbation theory. Weshow that the effects of transients on the kurtosis from the initial conditions, based onsecond-order Lagrangian perturbation theory (2LPT), have almost disappeared byz∼5, as long as the initialconditions are set at z>30. This means that for practical purposes, the initial conditions based on 2LPT are accurateenough for numerical calculations of skewness and kurtosis.

Initial Conditions to CosmologicalN‐Body Simulations, or, How to Run an Ensemble of Simulations

The conventional method of generating initial conditions for cosmological N-body simulations introduces a significant error in the real-space statistical properties of the particles. More specifically, the finite box size leads to a significant underestimate of σ8, the correlation function, and nonlinear effects. I implement a method of generating initial conditions for N-body simulations that accurately models the real-space statistical properties, such as the mass variance in spheres and the correlation function. The method requires running ensembles of simulations because the power in the DC mode is no longer assumed to be zero. For moderately sized boxes, I demonstrate that the new method corrects the underestimate in the mass variance in spheres and the shape of the correlation function. I also argue that subtracting Poisson noise from the power spectrum is a dangerous practice. Code to generate initial conditions to second order in Lagrangian perturbation theory with the new method is available from the author.

A Comprehensive Study of Young Black Hole Populations

We present theoretical models of black hole (BH) populations in young stellar environments, such as starbursts and young star clusters. Using a population synthesis approach, we compute the formation rates and characteristic properties of single and binary BHs for various representative ages and choices of parameters. We find that most of the BHs (typically 80% for an initial 50% binary fraction) are single, but with many originating from primordial binaries (which either merged into a single massive star or were disrupted following a supernova explosion). A smaller but significant fraction (typically 20%) of the BHs remain in binary systems. Main-sequence stars are the most frequent BH companions, but massive BH-BH binaries are the next most numerous group. The most massive BHs found in our simulations reach ~80 M☉ and are formed through mergers of massive binary components. If formed in a dense star cluster, such a massive stellar BH may become the seed for growth to a more massive (intermediate-mass) BH. Although we do not include dynamical interactions, our results provide realistic initial conditions for N-body simulations of dense star clusters (e.g., globular clusters) including primordial BHs.