SPH With Inter-dependent Fine-grained Tasking

Publications describing SWIFT

SWIFT: A modern highly-parallel gravity and smoothed particle hydrodynamics solver for astrophysical and cosmological applications

  Schaller M. et al., arXiv, 2023, 2305.13380   (citations: 12)

SWIFT: Using Task-Based Parallelism, Fully Asynchronous Communication, and Graph Partition-Based Domain Decomposition for Strong Scaling on more than 100,000 Cores

  Schaller M. et al., PASC, 2016, vol. 1

Publications using SWIFT

The following 58 publications have used SWIFT (directly or indirectly) to obtain their results. They have jointly gathered 560 citations.

58) The nature of compact radio sources: the case of FR 0 radio galaxies

  Baldi, R, A&ARv, 2023 , vol. 31 , issue 1   (citations: 0)


Radio-loud compact radio sources (CRSs) are characterised by morphological compactness of the jet structure centred on the active nucleus of the galaxy. Most of the local elliptical galaxies are found to host a CRS with nuclear luminosities lower than those of typical quasars, ≲1042ergs-1 . Recently, low-luminosity CRSs with a LINER-like optical spectrum have been named Fanaroff-Riley (FR) type 0 to highlight their lack of substantially extended radio emission at kpc scales, in contrast with the other Fanaroff-Riley classes, full-fledged FR Is and FR II radio galaxies. FR 0s are the most abundant class of radio galaxies in the local Universe, and characterised by a higher core dominance, poorer Mpc-scale environment and smaller (sub-kpc scale, if resolved) jets than FR Is. However, FR 0s share similar host and nuclear properties with FR Is. A different accretion-ejection paradigm from that in place in FR Is is invoked to account for the parsec-scale FR 0 jets. This review revises the state-of-the-art knowledge about FR 0s, their nature, and which open issues the next generation of radio telescopes can solve in this context.

57) Fisher matrix forecasts on the astrophysics of galaxies during the epoch of reionization from the 21-cm power spectra

  Balu, S et al., MNRAS, 2023 , vol. 525 , issue 2   (citations: 0)


The hyperfine 21-cm transition of neutral hydrogen from the early Universe (z > 5) is a sensitive probe of the formation and evolution of the first luminous sources. Using the Fisher matrix formalism we explore the complex and degenerate high-dimensional parameter space associated with the high-z sources of this era and forecast quantitative constraints from a future 21-cm power spectrum (21-cm PS) detection. This is achieved using $\rm {\small ERAXES}$, a coupled semi-analytic galaxy formation model and reionization simulation, applied to an N-body halo merger tree with a statistically complete population of all atomically cooled galaxies out to z ~ 20. Our mock observation assumes a 21-cm detection spanning z ∈ [5, 24] from a 1000 h mock observation with the forthcoming Square Kilometre Array, and is calibrated with respect to ultraviolet luminosity functions (UV LFs) at z ∈ [5, 10], the optical depth of CMB photons to Thompson scattering from Planck, and various constraints on the IGM neutral fraction at z > 5. In this work, we focus on the X-ray luminosity, ionizing UV photon escape fraction, star formation, and supernova feedback of the first galaxies. We demonstrate that it is possible to recover five of the eight parameters describing these properties with better than 50 per cent precision using just the 21-cm PS. By combining with UV LFs, we are able to improve our forecast, with five of the eight parameters constrained to better than 10 per cent (and all below 50 per cent).

56) A recent impact origin of Saturn's rings and mid-sized moons

  Teodoro, L et al., arXiv, 2023   (citations: 0)


We simulate the collision of precursor icy moons analogous to Dione and Rhea as a possible origin for Saturn's remarkably young rings. Such an event could have been triggered a few hundred million years ago by resonant instabilities in a previous satellite system. Using high-resolution smoothed particle hydrodynamics simulations, we find that this kind of impact can produce a wide distribution of massive objects and scatter material throughout the system. This includes the direct placement of pure-ice ejecta onto orbits that enter Saturn's Roche limit, which could form or rejuvenate rings. In addition, fragments and debris of rock and ice totalling more than the mass of Enceladus can be placed onto highly eccentric orbits that would intersect with any precursor moons orbiting in the vicinity of Mimas, Enceladus, or Tethys. This could prompt further disruption and facilitate a collisional cascade to distribute more debris for potential ring formation, the re-formation of the present-day moons, and evolution into an eventual cratering population of planeto-centric impactors.

55) Atmospheric loss in giant impacts depends on pre-impact surface conditions

  Lock, S et al., arXiv, 2023   (citations: 0)


Earth likely acquired much of its inventory of volatile elements during the main stage of its formation. Some of Earth's proto-atmosphere must therefore have survived the giant impacts, collisions between planet-sized bodies, that dominate the latter phases of accretion. Here we use a suite of 1D hydrodynamic simulations and impedance match calculations to quantify the effect that pre-impact surface conditions (such as atmospheric pressure and presence of an ocean) have on the efficiency of atmospheric and ocean loss from proto-planets during giant impacts. We find that -- in the absence of an ocean -- lighter, hotter, and lower-pressure atmospheres are more easily lost. The presence of an ocean can significantly increase the efficiency of atmospheric loss compared to the no-ocean case, with a rapid transition between low and high loss regimes as the mass ratio of atmosphere to ocean decreases. However, contrary to previous thinking, the presence of an ocean can also reduce atmospheric loss if the ocean is not sufficiently massive, typically less than a few times the atmospheric mass. Volatile loss due to giant impacts is thus highly sensitive to the surface conditions on the colliding bodies. To allow our results to be combined with 3D impact simulations, we have developed scaling laws that relate loss to the ground velocity and surface conditions. Our results demonstrate that the final volatile budgets of planets are critically dependent on the exact timing and sequence of impacts experienced by their precursor planetary embryos, making atmospheric properties a highly stochastic outcome of accretion.

54) First light and reionisation epoch simulations (FLARES) IX: The physical mechanisms driving compact galaxy formation and evolution

  Roper, W et al., MNRAS, 2023   (citations: 5)


In the FLARES (First Light And Reionisation Epoch Simulations) suite of hydrodynamical simulations, we find the high redshift (z > 5) intrinsic size-luminosity relation is, surprisingly, negatively sloped. However, after including the effects of dust attenuation we find a positively sloped UV observed size-luminosity relation in good agreement with other simulated and observational studies. In this work, we extend this analysis to probe the underlying physical mechanisms driving the formation and evolution of the compact galaxies driving the negative size-mass/size-luminosity relation. We find the majority of compact galaxies (R1/2, ⋆ < 1 physical kpc, which drive the negative slope of the size-mass relation, have transitioned from extended to compact sizes via efficient centralised cooling, resulting in high specific star formation rates in their cores. These compact stellar systems are enshrouded by non-star forming gas distributions as much as 100 × larger than their stellar counterparts. By comparing with galaxies from the EAGLE simulation suite, we find that these extended gas distributions 'turn on' and begin to form stars between z = 5 and z = 0 leading to increasing sizes, and thus the evolution of the size-mass relation from a negative to a positive slope. This explicitly demonstrates the process of inside-out galaxy formation in which compact bulges form earlier than the surrounding discs.

53) The FLAMINGO project: revisiting the $S_8$ tension and the role of baryonic physics

  McCarthy, I et al., arXiv, 2023   (citations: 0)


A number of recent studies have found evidence for a tension between observations of large-scale structure (LSS) and the predictions of the standard model of cosmology with the cosmological parameters fit to the cosmic microwave background (CMB). The origin of this '$S_8$ tension' remains unclear, but possibilities include new physics beyond the standard model, unaccounted for systematic errors in the observational measurements and/or uncertainties in the role that baryons play. Here we carefully examine the latter possibility using the new FLAMINGO suite of large-volume cosmological hydrodynamical simulations. We project the simulations onto observable harmonic space and compare with observational measurements of the power and cross-power spectra of cosmic shear, CMB lensing, and the thermal Sunyaev-Zel'dovich (tSZ) effect. We explore the dependence of the predictions on box size and resolution, cosmological parameters including the neutrino mass, and the efficiency and nature of baryonic 'feedback'. Despite the wide range of astrophysical behaviours simulated, we find that baryonic effects are not sufficiently large to remove the $S_8$ tension. Consistent with recent studies, we find the CMB lensing power spectrum is in excellent agreement with the standard model, whilst the cosmic shear power spectrum, tSZ effect power spectrum, and the cross-spectra between shear, CMB lensing, and the tSZ effect are all in varying degrees of tension with the CMB-specified standard model. These results suggest that some mechanism is required to slow the growth of fluctuations at late times and/or on non-linear scales, but that it is unlikely that baryon physics is driving this modification.

52) Tests of subgrid models for star formation using simulations of isolated disk galaxies

  Nobels, F et al., arXiv, 2023   (citations: 0)


We use smoothed-particle hydrodynamics simulations of isolated Milky Way-mass disk galaxies that include cold, interstellar gas to test subgrid prescriptions for star formation (SF). Our fiducial model combines a Schmidt law with a gravitational instability criterion, but we also test density thresholds and temperature ceilings. While SF histories are insensitive to the prescription for SF, the Kennicutt-Schmidt (KS) relations between SF rate and gas surface density can discriminate between models. We show that our fiducial model, with an SF efficiency per free-fall time of 1 per cent, agrees with spatially-resolved and azimuthally-averaged observed KS relations for neutral, atomic and molecular gas. Density thresholds do not perform as well. While temperature ceilings selecting cold, molecular gas can match the data for galaxies with solar metallicity, they are unsuitable for very low-metallicity gas and hence for cosmological simulations. We argue that SF criteria should be applied at the resolution limit rather than at a fixed physical scale, which means that we should aim for numerical convergence of observables rather than of the properties of gas labelled as star-forming. Our fiducial model yields good convergence when the mass resolution is varied by nearly 4 orders of magnitude, with the exception of the spatially-resolved molecular KS relation at low surface densities. For the gravitational instability criterion, we quantify the impact on the KS relations of gravitational softening, the SF efficiency, and the strength of supernova feedback, as well as of observable parameters such as the inclusion of ionized gas, the averaging scale, and the metallicity.

51) The FLAMINGO project: cosmological hydrodynamical simulations for large-scale structure and galaxy cluster surveys

  Schaye, J et al., MNRAS, 2023   (citations: 17)


We introduce the Virgo Consortium's FLAMINGO suite of hydrodynamical simulations for cosmology and galaxy cluster physics. To ensure the simulations are sufficiently realistic for studies of large-scale structure, the subgrid prescriptions for stellar and AGN feedback are calibrated to the observed low-redshift galaxy stellar mass function and cluster gas fractions. The calibration is performed using machine learning, separately for each of FLAMINGO's three resolutions. This approach enables specification of the model by the observables to which they are calibrated. The calibration accounts for a number of potential observational biases and for random errors in the observed stellar masses. The two most demanding simulations have box sizes of 1.0 and 2.8 Gpc on a side and baryonic particle masses of 1 × 108 and 1 × 109 M, respectively. For the latter resolution the suite includes 12 model variations in a 1 Gpc box. There are 8 variations at fixed cosmology, including shifts in the stellar mass function and/or the cluster gas fractions to which we calibrate, and two alternative implementations of AGN feedback (thermal or jets). The remaining 4 variations use the unmodified calibration data but different cosmologies, including different neutrino masses. The 2.8 Gpc simulation follows 3 × 1011 particles, making it the largest ever hydrodynamical simulation run to z = 0. Lightcone output is produced on-the-fly for up to 8 different observers. We investigate numerical convergence, show that the simulations reproduce the calibration data, and compare with a number of galaxy, cluster, and large-scale structure observations, finding very good agreement with the data for converged predictions. Finally, by comparing hydrodynamical and 'dark-matter-only' simulations, we confirm that baryonic effects can suppress the halo mass function and the matter power spectrum by up to ≈20 per cent.

50) A thermal-kinetic subgrid model for supernova feedback in simulations of galaxy formation

  Chaikin, E et al., MNRAS, 2023 , vol. 523 , issue 3   (citations: 5)


We present a subgrid model for supernova feedback designed for cosmological simulations of galaxy formation that may include a cold interstellar medium (ISM). The model uses thermal and kinetic channels of energy injection, which are built upon the stochastic kinetic and thermal models for stellar feedback used in the OWLS and EAGLE simulations, respectively. In the thermal channel, the energy is distributed statistically isotropically and injected stochastically in large amounts per event, which minimizes spurious radiative energy losses. In the kinetic channel, we inject the energy in small portions by kicking gas particles in pairs in opposite directions. The implementation of kinetic feedback is designed to conserve energy, linear and angular momentum, and is statistically isotropic. To test the model, we run simulations of isolated Milky Way-mass and dwarf galaxies, in which the gas is allowed to cool down to 10 K. Using the thermal and kinetic channels together, we obtain smooth star formation histories and powerful galactic winds with realistic mass loading factors. Furthermore, the model produces spatially resolved star formation rates (SFRs) and velocity dispersions that are in agreement with observations. We vary the numerical resolution by several orders of magnitude and find excellent convergence of the global SFRs and wind mass loading. We show that large thermal energy injections generate a hot phase of the ISM and modulate the star formation by ejecting gas from the disc, while the low-energy kicks increase the turbulent velocity dispersion in the neutral ISM, which in turn helps suppress star formation.

49) Implications of z ≳ 12 JWST galaxies for galaxy formation at high redshift

  Qin, Y et al., MNRAS, 2023   (citations: 5)


Using a semi-analytic galaxy-formation model, we study analogues of 8 recently discovered JWST galaxies at z ≳ 12. We select analogues from a cosmological simulation with a (311cMpc)3 volume and an effective particle number of 1012 enabling resolution of every atomic-cooling galaxy at z ≤ 20. We vary model parameters to reproduce the observed UV luminosity function at 5 < z < 13, aiming for a statistically representative high-redshift galaxy mock catalogue. Using the forward-modelled JWST photometry, we identify analogues from this catalogue and study their properties as well as possible evolutionary paths and local environment. We find faint JWST galaxies (MUV ≳ - 19.5) to remain consistent with standard galaxy-formation model and that our fiducial catalogue includes large samples of their analogues. The properties of these analogues broadly agree with conventional SED fitting results, except for having systematically lower redshifts due to the evolving UV luminosity function, and for having higher specific star formation rates as a result of burstier histories in our model. On the other hand, only a handful of bright galaxy analogues can be identified for the observed z ~ 12 galaxies. Moreover, in order to reproduce the z ≳ 16 JWST galaxy candidates, boosting star-forming efficiencies through reduced feedback regulation and increased gas depletion rate is necessary relative to models of lower-redshift populations. This suggests star formation in the first galaxies could differ significantly from their lower-redshift counterparts. We also find that these candidates are subject to low-redshift contamination, which is present in our fiducial results as both the dusty or quiescent galaxies at z ~ 5.

48) The compactness of ultra faint dwarf galaxies : a new challenge ?

  Revaz, Y, arXiv, 2023   (citations: 0)


So far, numerical simulations of ultra-faint dwarf galaxies (UFDs) failed to properly reproduce the observed size-luminosity relation. In particular, no hydro-dynamical-run managed to form UFDs with a half light radius as small as 30 pc as seen in several UFD candidates. We tackle this problem by developing a simple but numerically clean and powerful method in which predictions of the stellar content of UFDs from LCDM cosmological hydro-dynamical-simulations is combined with very high resolution dark matter only runs. This method allows to trace the build-up history of UFDs and determine the impact of the merger of building-block objects on their final size. We found that, while no UFDs more compact than 20 pc can be formed, slightly larger system are reproduced only if all member stars are issued from the same initial mini-halo. However this imposes (i) the total virial mass to be smaller than 3x10^8Msol, (ii) the stellar content prior to the end of the re-ionisation epoch to be very compact (<15 pc) and strongly gravitationally bound, a challenge for current hydro-dynamical numerical simulations. If initial stellar building blocks are larger than 35 pc the size of the UFD will extend to 80 pc. Finally, our study shows that UFDs keep strong imprints of their build-up history in the form of elongated or extended stellar halos. Those features can erroneously be interpreted as tidal signatures.

47) Multi-epoch machine learning 2: identifying physical drivers of galaxy properties in simulations

  McGibbon, R et al., MNRAS, 2023 , vol. 523 , issue 4   (citations: 0)


Using a novel machine learning method, we investigate the buildup of galaxy properties in different simulations, and in various environments within a single simulation. The aim of this work is to show the power of this approach at identifying the physical drivers of galaxy properties within simulations. We compare how the stellar mass is dependent on the value of other galaxy and halo properties at different points in time by examining the feature importance values of a machine learning model. By training the model on IllustrisTNG, we show that stars are produced at earlier times in higher density regions of the universe than they are in low density regions. We also apply the technique to the Illustris, EAGLE, and CAMELS simulations. We find that stellar mass is built up in a similar way in EAGLE and IllustrisTNG, but significantly differently in the original Illustris, suggesting that subgrid model physics is more important than the choice of hydrodynamics method. These differences are driven by the efficiency of supernova feedback. Applying principal component analysis to the CAMELS simulations allows us to identify a component associated with the importance of a halo's gravitational potential and another component representing the time at which galaxies form. We discover that the speed of galactic winds is a more critical subgrid parameter than the total energy per unit star formation. Finally, we find that the Simba black hole feedback model has a larger effect on galaxy formation than the IllustrisTNG black hole feedback model.

46) Sensitivity of non-radiative cloud-wind interactions to the hydrodynamic solver

  Braspenning, J et al., MNRAS, 2023 , vol. 523 , issue 1   (citations: 6)


Cloud-wind interactions are common in the interstellar and circumgalactic media. Many studies have used simulations of such interactions to investigate the effect of particular physical processes, but the impact of the choice of hydrodynamic solver has largely been overlooked. Here we study the cloud-wind interaction, also known as the 'blob test', using seven different hydrodynamic solvers: three flavours of SPH, a moving mesh, adaptive mesh refinement, and two meshless schemes. The evolution of masses in dense gas and intermediate-temperature gas, as well as the covering fraction of intermediate-temperature gas, are systematically compared for initial density contrasts of 10 and 100, and five numerical resolutions. To isolate the differences due to the hydrodynamic solvers, we use idealized non-radiative simulations without physical conduction. We find large differences between these methods. SPH methods show slower dispersal of the cloud, particularly for the higher density contrast, but faster convergence, especially for the lower density contrast. Predictions for the intermediate-temperature gas differ particularly strongly, also between non-SPH codes, and converge most slowly. We conclude that the hydrodynamical interaction between a dense cloud and a supersonic wind remains an unsolved problem. Studies aiming to understand the physics or observational signatures of cloud-wind interactions should test the robustness of their results by comparing different hydrodynamic solvers.

45) Robust Field-level Likelihood-free Inference with Galaxies

  de Santi, N et al., ApJ, 2023 , vol. 952 , issue 1   (citations: 8)


We train graph neural networks to perform field-level likelihood-free inference using galaxy catalogs from state-of-the-art hydrodynamic simulations of the CAMELS project. Our models are rotational, translational, and permutation invariant and do not impose any cut on scale. From galaxy catalogs that only contain 3D positions and radial velocities of ~1000 galaxies in tiny ${(25\,{h}^{-1}\mathrm{Mpc})}^{3}$ volumes our models can infer the value of Ωm with approximately 12% precision. More importantly, by testing the models on galaxy catalogs from thousands of hydrodynamic simulations, each having a different efficiency of supernova and active galactic nucleus feedback, run with five different codes and subgrid models-IllustrisTNG, SIMBA, Astrid, Magneticum, SWIFT-EAGLE-we find that our models are robust to changes in astrophysics, subgrid physics, and subhalo/galaxy finder. Furthermore, we test our models on 1024 simulations that cover a vast region in parameter space-variations in five cosmological and 23 astrophysical parameters-finding that the model extrapolates really well. Our results indicate that the key to building a robust model is the use of both galaxy positions and velocities, suggesting that the network has likely learned an underlying physical relation that does not depend on galaxy formation and is valid on scales larger than ~10 h -1 kpc.

44) Winds versus jets: a comparison between black hole feedback modes in simulations of idealized galaxy groups and clusters

  Husko, F et al., arXiv, 2023   (citations: 1)


Using the SWIFT simulation code we study different forms of active galactic nuclei (AGN) feedback in idealized galaxy groups and clusters. We first present a physically motivated model of black hole (BH) spin evolution and a numerical implementation of thermal isotropic feedback (representing the effects of energy-driven winds) and collimated kinetic jets that they launch at different accretion rates. We find that kinetic jet feedback is more efficient at quenching star formation in the brightest cluster galaxies (BCGs) than thermal isotropic feedback, while simultaneously yielding cooler cores in the intracluster medium (ICM). A hybrid model with both types of AGN feedback yields moderate star formation rates, while having the coolest cores. We then consider a simplified implementation of AGN feedback by fixing the feedback efficiencies and the jet direction, finding that the same general conclusions hold. We vary the feedback energetics (the kick velocity and the heating temperature), the fixed efficiencies and the type of energy (kinetic versus thermal) in both the isotropic and the jet case. The isotropic case is largely insensitive to these variations. In particular, we highlight that kinetic isotropic feedback (used e.g. in IllustrisTNG) is similar in its effects to its thermal counterpart (used e.g. in EAGLE). On the other hand, jet feedback must be kinetic in order to be efficient at quenching. We also find that it is much more sensitive to the choice of energy per feedback event (the jet velocity), as well as the efficiency. The former indicates that jet velocities need to be carefully chosen in cosmological simulations, while the latter motivates the use of BH spin evolution models.

43) Influence of local structure on relic neutrino abundances and anisotropies

  Zimmer, F et al., arXiv, 2023   (citations: 2)


Gravitational potentials of the Milky Way and extragalactic structures can influence the propagation of the cosmic neutrino background (CNB). Of particular interest to future CNB observatories, such as PTOLEMY, is the CNB number density on Earth. In this study, we have developed a simulation framework that maps the trajectories of relic neutrinos as they move through the local gravitational environment. The potentials are based on the dark matter halos found in state-of-the-art cosmological N-body simulations, resulting in a more nuanced and realistic input than the previously employed analytical models. We find that the complex dark matter distributions, along with their dynamic evolution, influence the abundance and anisotropies of the CNB in ways unaccounted for by earlier analytical methods. Importantly, these cosmological simulations contain multiple instances of Milky Way-like halos that we employ to model a variety of gravitational landscapes. Consequently, we notice a variation in the CNB number densities that can be primarily attributed to the differences in the masses of these individual halos. For neutrino masses between $0.01$ and $0.3$ eV, we note clustering factors within the range of $1+\mathcal{O}(10^{-3})$ to $1+\mathcal{O}(1)$. Furthermore, the asymmetric nature of the underlying dark matter distributions within the halos results in not only overdense, but intriguingly, underdense regions within the full-sky anisotropy maps. Gravitational clustering appears to have a significant impact on the angular power spectra of these maps, leading to orders of magnitude more power on smaller scales beyond multipoles of $\ell = 3$ when juxtaposed against predictions by primordial fluctuations. We discuss how our results reshape our understanding of relic neutrino clustering and how this might affect observability of future CNB observatories such as PTOLEMY.

42) Sarracen: a Python package for analysis and visualization of smoothed particle hydrodynamics data

  Harris, A et al., JOSS, 2023 , vol. 8 , issue 86   (citations: 1)

    No abstract

41) FLAMINGO: Calibrating large cosmological hydrodynamical simulations with machine learning

  Kugel, R et al., arXiv, 2023   (citations: 4)


To fully take advantage of the data provided by large-scale structure surveys, we need to quantify the potential impact of baryonic effects, such as feedback from active galactic nuclei (AGN) and star formation, on cosmological observables. In simulations, feedback processes originate on scales that remain unresolved. Therefore, they need to be sourced via subgrid models that contain free parameters. We use machine learning to calibrate the AGN and stellar feedback models for the FLAMINGO cosmological hydrodynamical simulations. Using Gaussian process emulators trained on Latin hypercubes of 32 smaller-volume simulations, we model how the galaxy stellar mass function and cluster gas fractions change as a function of the subgrid parameters. The emulators are then fit to observational data, allowing for the inclusion of potential observational biases. We apply our method to the three different FLAMINGO resolutions, spanning a factor of 64 in particle mass, recovering the observed relations within the respective resolved mass ranges. We also use the emulators, which link changes in subgrid parameters to changes in observables, to find models that skirt or exceed the observationally allowed range for cluster gas fractions and the stellar mass function. Our method enables us to define model variations in terms of the data that they are calibrated to rather than the values of specific subgrid parameters. This approach is useful, because subgrid parameters are typically not directly linked to particular observables, and predictions for a specific observable are influenced by multiple subgrid parameters.

40) Euclid: modelling massive neutrinos in cosmology - a code comparison

  Adamek, J et al., JCAP, 2023 , vol. 2023 , issue 6   (citations: 11)


The measurement of the absolute neutrino mass scale from cosmological large-scale clustering data is one of the key science goals of the Euclid mission. Such a measurement relies on precise modelling of the impact of neutrinos on structure formation, which can be studied with N-body simulations. Here we present the results from a major code comparison effort to establish the maturity and reliability of numerical methods for treating massive neutrinos. The comparison includes eleven full N-body implementations (not all of them independent), two N-body schemes with approximate time integration, and four additional codes that directly predict or emulate the matter power spectrum. Using a common set of initial data we quantify the relative agreement on the nonlinear power spectrum of cold dark matter and baryons and, for the N-body codes, also the relative agreement on the bispectrum, halo mass function, and halo bias. We find that the different numerical implementations produce fully consistent results. We can therefore be confident that we can model the impact of massive neutrinos at the sub-percent level in the most common summary statistics. We also provide a code validation pipeline for future reference.

39) Testing Jeans dynamical models with prolate rotation on a cosmologically simulated dwarf galaxy

  Sedain, A et al., arXiv, 2023   (citations: 0)


Prolate rotation is characterized by a significant stellar rotation around a galaxy's major axis, which contrasts with the more common oblate rotation. Prolate rotation is thought to be due to major mergers and thus studies of prolate-rotating systems can help us better understand the hierarchical process of galaxy evolution. Dynamical studies of such galaxies are important to find their gravitational potential profile, total mass, and dark matter fraction. Recently, it has been shown in a cosmological simulation that it is possible to form a prolate-rotating dwarf galaxy following a dwarf-dwarf merger event. The simulation also shows that the unusual prolate rotation can be time enduring. In this particular example, the galaxy continued to rotate around its major axis for at least $7.4$\,Gyr (from the merger event until the end of the simulation). In this project, we use mock observations of the hydro-dynamically simulated prolate-rotating dwarf galaxy to fit various stages of its evolution with Jeans dynamical models. The Jeans models successfully fit the early oblate state before the major merger event, and also the late prolate stages of the simulated galaxy, recovering its mass distribution, velocity dispersion, and rotation profile. We also ran a prolate-rotating N-body simulation with similar properties to the cosmologically simulated galaxy, which gradually loses its angular momentum on a short time scale $\sim100$\,Myr. More tests are needed to understand why prolate rotation is time enduring in the cosmological simulation, but not in a simple N-body simulation.

38) The complex interplay of AGN jet-inflated bubbles and the intracluster medium

  Husko, F et al., MNRAS, 2023 , vol. 521 , issue 3   (citations: 6)


We use SWIFT, a smoothed particle hydrodynamics code, to simulate the evolution of bubbles inflated by active galactic nuclei (AGNs) jets, as well as their interactions with the ambient intracluster medium (ICM). These jets inflate lobes that turn into bubbles after the jets are turned off (at t = 50 Myr). Almost all of the energy injected into the jets is transferred to the ICM very quickly after they are turned off, with roughly 70 per cent of it in thermal form and the rest in kinetic. At late times (t > 500 Myr) we find the following: (1) the bubbles draw out trailing filaments of low-entropy gas, similar to those recently observed, (2) the action of buoyancy and the uplift of the filaments dominates the energetics of both the bubbles and the ICM, and (3) almost all of the originally injected energy is in the form of gravitational potential energy, with the bubbles containing 15 per cent of it, and the rest contained in the ICM. These findings indicate that feedback proceeds mainly through the displacement of gas to larger radii. We find that the uplift of these filaments permanently changes the thermodynamic properties of the ICM by reducing the central density and increasing the central temperature (within 30 kpc). We propose that jet feedback proceeds not only through the heating of the ICM (which can delay cooling), but also through the uplift-related reduction of the central gas density. The latter also delays cooling, on top of reducing the amount of gas available to cool.

37) The impact and response of minihalos and the inter-halo medium on cosmic reionization

  Chan, T et al., arXiv, 2023   (citations: 1)


An ionization front (I-front) that propagates through an inhomogeneous medium is slowed down by self-shielding and recombinations. We perform cosmological radiation hydrodynamics simulations of the I-front propagation during the epoch of cosmic reionization. The simulations resolve gas in minihalos (halo mass $10^4\lesssim M_h[{\rm M}_\odot]\lesssim 10^8)$ that could dominate recombinations, in a computational volume that is large enough to sample the abundance of such halos. The numerical resolution is sufficient (gas particle mass $\sim 20{\rm M}_\odot$, spatial resolution $< 0.1\;{\rm ckpc}$) to allow accurate modelling of the hydrodynamic response of gas to photo-heating. We quantify the photo-evaporation time of minihalos as a function of $M_h$ and its dependence on the photo-ionization rate, $\Gamma_{-12}$, and the redshift of reionization, $z_i$. The recombination rate can be enhanced over that of a uniform medium by a factor $\sim 10-20$ early on. The peak value increases with $\Gamma_{-12}$ and decreases with $z_i$, due to the enhanced contribution from minihalos. The clumping factor, $c_r$, decreases to a factor of a few at $\sim 100\;{\rm Myr}$ after the passage of the I-front when the minihalos have been photo-evaporated; this asymptotic value depends only weakly on $\Gamma_{-12}$. Recombinations increase the required number of photons per baryon to reionize the Universe by 20-100 per cent, with the higher value occurring when $\Gamma_{-12}$ is high and $z_i$ is low. We complement the numerical simulations with simple analytical models for the evaporation rate and the inverse Strömgren layer. The study also demonstrates the proficiency and potential of SPHM1RT to address astrophysical problems in high-resolution cosmological simulations.

36) EAGLE-like simulation models do not solve the entropy core problem in groups and clusters of galaxies

  Altamura, E et al., MNRAS, 2023 , vol. 520 , issue 2   (citations: 4)


Recent high-resolution cosmological hydrodynamic simulations run with a variety of codes systematically predict large amounts of entropy in the intra-cluster medium at low redshift, leading to flat entropy profiles and a suppressed cool-core population. This prediction is at odds with X-ray observations of groups and clusters. We use a new implementation of the EAGLE galaxy formation model to investigate the sensitivity of the central entropy and the shape of the profiles to changes in the sub-grid model applied to a suite of zoom-in cosmological simulations of a group of mass M500 = 8.8 × 1012 M and a cluster of mass 2.9 × 1014 M. Using our reference model, calibrated to match the stellar mass function of field galaxies, we confirm that our simulated groups and clusters contain hot gas with too high entropy in their cores. Additional simulations run without artificial conduction, metal cooling or active galactic nuclei (AGN) feedback produce lower entropy levels but still fail to reproduce observed profiles. Conversely, the two objects run without supernova feedback show a significant entropy increase which can be attributed to excessive cooling and star formation. Varying the AGN heating temperature does not greatly affect the profile shape, but only the overall normalization. Finally, we compared runs with four AGN heating schemes and obtained similar profiles, with the exception of bipolar AGN heating, which produces a higher and more uniform entropy distribution. Our study leaves open the question of whether the entropy core problem in simulations, and particularly the lack of power-law cool-core profiles, arise from incorrect physical assumptions, missing physical processes, or insufficient numerical resolution.

35) Active galactic nuclei jets simulated with smoothed particle hydrodynamics

  Husko, F et al., MNRAS, 2023 , vol. 520 , issue 4   (citations: 8)


Simulations of active galactic nuclei (AGN) jets have thus far been performed almost exclusively using grid-based codes. We present the first results from hydrodynamical tests of AGN jets, and their interaction with the intracluster medium (ICM), using smoothed particle hydrodynamics as implemented in the SWIFT code. We launch these jets into a constant-density ICM, as well as ones with a power-law density profile. We also vary the jet power, velocity, opening angle, and numerical resolution. In all cases we find broad agreement between our jets and theoretical predictions for the lengths of the jets and the lobes they inflate, as well as the radii of the lobes. The jets first evolve ballistically, and then transition to a self-similar phase, during which the lobes expand in a self-similar fashion (keeping a constant shape). In this phase the kinetic and thermal energies in the lobes and in the shocked ICM are constant fractions of the total injected energy. In our standard simulation, two thirds of the initially injected energy is transferred to the ICM by the time the jets are turned off, mainly through a bow shock. Of that, $70{{\%}}$ is in kinetic form, indicating that the bow shock does not fully and efficiently thermalize while the jet is active. At resolutions typical of large cosmological simulations (mgas ≈ 107 M), the shape of the lobes is close to self-similar predictions to an accuracy of $15{{\%}}$. This indicates that the basic physics of jet-inflated lobes can be correctly simulated even at such resolutions (≈500 particles per jet).

34) Galaxy Evolution in $\ddot{\mu}$ based Cosmologies

  Roper, W et al., arXiv, 2023   (citations: 0)


We present the first study of galaxy evolution in $\ddot{\mu}$ based cosmologies. We find that recent JWST observations of massive galaxies at extremely high redshifts are consistent with such a cosmology. However, the low redshift Universe is entirely divergent from the $\ddot{\mu}$ cosmic star formation rate density. We thus propose that our Universe was at one point dominated by a Primordial Bovine Herd (PBH) which later decayed producing dark energy. Note that we do not detail the mechanisms by which this decay process takes place. Despite its vanishingly small probability for existence, a $\ddot{\mu}$ based cosmological model marries the disparate findings in the high and low redshift Universe.

33) Mesh-free hydrodynamics in PKDGRAV3 for galaxy formation simulations

  Alonso Asensio, I et al., MNRAS, 2023 , vol. 519 , issue 1   (citations: 4)


We extend the state-of-the-art N-body code PKDGRAV3 with the inclusion of mesh-free gas hydrodynamics for cosmological simulations. Two new hydrodynamic solvers have been implemented, the mesh-less finite volume and mesh-less finite mass methods. The solvers manifestly conserve mass, momentum, and energy, and have been validated with a wide range of standard test simulations, including cosmological simulations. We also describe improvements to PKDGRAV3 that have been implemented for performing hydrodynamic simulations. These changes have been made with efficiency and modularity in mind, and provide a solid base for the implementation of the required modules for galaxy formation and evolution physics and future porting to GPUs. The code is released in a public repository, together with the documentation, and all the test simulations presented in this work.

32) GEAR-RT: Towards Exa-Scale Moment Based Radiative Transfer For Cosmological Simulations Using Task-Based Parallelism And Dynamic Sub-Cycling with SWIFT

  Ivkovic, M, arXiv, 2023   (citations: 2)


The development and implementation of GEAR-RT, a radiative transfer solver using the M1 closure in the open source code SWIFT, is presented, and validated using standard tests for radiative transfer. GEAR-RT is modeled after RAMSES-RT (Rosdahl et al. 2013) with some key differences. Firstly, while RAMSES-RT uses Finite Volume methods and an Adaptive Mesh Refinement (AMR) strategy, GEAR-RT employs particles as discretization elements and solves the equations using a Finite Volume Particle Method (FVPM). Secondly, GEAR-RT makes use of the task-based parallelization strategy of SWIFT, which allows for optimized load balancing, increased cache efficiency, asynchronous communications, and a domain decomposition based on work rather than on data. GEAR-RT is able to perform sub-cycles of radiative transfer steps w.r.t. a single hydrodynamics step. Radiation requires much smaller time step sizes than hydrodynamics, and sub-cycling permits calculations which are not strictly necessary to be skipped. Indeed, in a test case with gravity, hydrodynamics, and radiative transfer, the sub-cycling is able to reduce the runtime of a simulation by over 90%. Allowing only a part of the involved physics to be sub-cycled is a contrived matter when task-based parallelism is involved, and is an entirely novel feature in SWIFT. Since GEAR-RT uses a FVPM, a detailed introduction into Finite Volume methods and Finite Volume Particle Methods is presented. In astrophysical literature, two FVPM methods are written about: Hopkins (2015) have implemented one in their GIZMO code, while the one mentioned in Ivanova et al. (2013) isn't used to date. In this work, I test an implementation of the Ivanova et al. (2013) version, and conclude that in its current form, it is not suitable for use with particles which are co-moving with the fluid, which in turn is an essential feature for cosmological simulations.

31) A universal equation to predict $\Omega_{\rm m}$ from halo and galaxy catalogues

  Shao, H et al., arXiv, 2023   (citations: 7)


We discover analytic equations that can infer the value of $\Omega_{\rm m}$ from the positions and velocity moduli of halo and galaxy catalogues. The equations are derived by combining a tailored graph neural network (GNN) architecture with symbolic regression. We first train the GNN on dark matter halos from Gadget N-body simulations to perform field-level likelihood-free inference, and show that our model can infer $\Omega_{\rm m}$ with $\sim6\%$ accuracy from halo catalogues of thousands of N-body simulations run with six different codes: Abacus, CUBEP$^3$M, Gadget, Enzo, PKDGrav3, and Ramses. By applying symbolic regression to the different parts comprising the GNN, we derive equations that can predict $\Omega_{\rm m}$ from halo catalogues of simulations run with all of the above codes with accuracies similar to those of the GNN. We show that by tuning a single free parameter, our equations can also infer the value of $\Omega_{\rm m}$ from galaxy catalogues of thousands of state-of-the-art hydrodynamic simulations of the CAMELS project, each with a different astrophysics model, run with five distinct codes that employ different subgrid physics: IllustrisTNG, SIMBA, Astrid, Magneticum, SWIFT-EAGLE. Furthermore, the equations also perform well when tested on galaxy catalogues from simulations covering a vast region in parameter space that samples variations in 5 cosmological and 23 astrophysical parameters. We speculate that the equations may reflect the existence of a fundamental physics relation between the phase-space distribution of generic tracers and $\Omega_{\rm m}$, one that is not affected by galaxy formation physics down to scales as small as $10~h^{-1}{\rm kpc}$.

30) Simulations of the reionization of the clumpy intergalactic medium with a novel particle-based two-moment radiative transfer scheme

  Chan, T et al., IAUS, 2023 , vol. 362   (citations: 1)


The progress of cosmic reionization depends on the presence of over-dense regions that act as photon sinks. Such sinks may slow down ionization fronts as compared to a uniform intergalactic medium (IGM) by increasing the clumping factor. We present simulations of reionization in a clumpy IGM resolving even the smallest sinks. The simulations use a novel, spatially adaptive and efficient radiative transfer implementation in the SWIFT SPH code, based on the two-moment method. We find that photon sinks can increase the clumping factor by a factor of ∼10 during the first ∼100 Myrs after the passage of an ionization front. After this time, the clumping factor decreases as the smaller sinks photoevaporate. Altogether, photon sinks increase the number of photons required to reionize the Universe by a factor of η ∼2, as compared to the homogeneous case. The value of η also depends on the emissivity of the ionizing sources.

29) TangoSIDM: tantalizing models of self-interacting dark matter

  Correa, C et al., MNRAS, 2022 , vol. 517 , issue 2   (citations: 11)


We introduce the TangoSIDM project, a suite of cosmological simulations of structure formation in a Λ-self-interacting dark matter (SIDM) universe. TangoSIDM explores the impact of large dark matter (DM) scattering cross-sections over dwarf galaxy scales. Motivated by DM interactions that follow a Yukawa potential, the cross-section per unit mass, σ/mχ, assumes a velocity-dependent form that avoids violations of current constraints on large scales. We demonstrate that our implementation accurately models not only core formation in haloes but also gravothermal core collapse. For central haloes in cosmological volumes, frequent DM particle collisions isotropise the particles orbit, making them largely spherical. We show that the velocity-dependent σ/mχ models produce a large diversity in the circular velocities of satellites haloes, with the spread in velocities increasing as the cross-sections reach 20, 60, and 100 cm2 g-1 in $10^9~\rm {M}_{\odot }$ haloes. The large variation in the haloes internal structure is driven by DM particles interactions, causing in some haloes the formation of extended cores, whereas in others gravothermal core collapse. We conclude that the SIDM models from the Tango project offer a promising explanation for the diversity in the density and velocity profiles of observed dwarf galaxies.

28) Higher order initial conditions with massive neutrinos

  Elbers, W et al., MNRAS, 2022 , vol. 516 , issue 3   (citations: 11)


The discovery that neutrinos have mass has important consequences for cosmology. The main effect of massive neutrinos is to suppress the growth of cosmic structure on small scales. Such growth can be accurately modelled using cosmological N-body simulations, but doing so requires accurate initial conditions (ICs). There is a trade-off, especially with first-order ICs, between truncation errors for late starts and discreteness and relativistic errors for early starts. Errors can be minimized by starting simulations at late times using higher order ICs. In this paper, we show that neutrino effects can be absorbed into scale-independent coefficients in higher order Lagrangian perturbation theory (LPT). This clears the way for the use of higher order ICs for massive neutrino simulations. We demonstrate that going to higher order substantially improves the accuracy of simulations. To match the sensitivity of surveys like DESI and Euclid, errors in the matter power spectrum should be well below $1{{\ \rm per\ cent}}$. However, we find that first-order Zel'dovich ICs lead to much larger errors, even when starting as early as z = 127, exceeding $1{{\ \rm per\ cent}}$ at z = 0 for k > 0.5 Mpc-1 for the power spectrum and k > 0.1 Mpc-1 for the equilateral bispectrum in our simulations. Ratios of power spectra with different neutrino masses are more robust than absolute statistics, but still depend on the choice of ICs. For all statistics considered, we obtain $1{{\ \rm per\ cent}}$ agreement between 2LPT and 3LPT at z = 0.

27) The impact of stochastic modeling on the predictive power of galaxy formation simulations

  Borrow, J et al., arXiv, 2022   (citations: 11)


All modern galaxy formation models employ stochastic elements in their sub-grid prescriptions to discretise continuous equations across the time domain. In this paper, we investigate how the stochastic nature of these models, notably star formation, black hole accretion, and their associated feedback, that act on small ($<$ kpc) scales, can back-react on macroscopic galaxy properties (e.g. stellar mass and size) across long ($>$ Gyr) timescales. We find that the scatter in scaling relations predicted by the EAGLE model implemented in the SWIFT code can be significantly impacted by random variability between re-simulations of the same object, even when galaxies are resolved by tens of thousands of particles. We then illustrate how re-simulations of the same object can be used to better understand the underlying model, by showing how correlations between galaxy stellar mass and black hole mass disappear at the highest black hole masses ($M_{\rm BH} > 10^8$ M$_\odot$), indicating that the feedback cycle may be interrupted by external processes. We find that although properties that are collected cumulatively over many objects are relatively robust against random variability (e.g. the median of a scaling relation), the properties of individual galaxies (such as galaxy stellar mass) can vary by up to 25\%, even far into the well-resolved regime, driven by bursty physics (black hole feedback) and mergers between galaxies. We suggest that studies of individual objects within cosmological simulations be treated with caution, and that any studies aiming to closely investigate such objects must account for random variability within their results.

26) Spin-driven jet feedback in idealized simulations of galaxy groups and clusters

  Husko, F et al., MNRAS, 2022 , vol. 516 , issue 3   (citations: 13)


We implement a black hole spin evolution and jet feedback model into SWIFT, a smoothed particle hydrodynamics code. The jet power is determined self-consistently assuming that the black hole accretion rate is equal to the Bondi rate (i.e. the accretion efficiency is 100 per cent), and using a realistic, spin-dependent efficiency. The jets are launched along the spin axis of the black hole, resulting in natural reorientation and precession. We apply the model to idealized simulations of galaxy groups and clusters, finding that jet feedback successfully quenches gas cooling and star formation in all systems. Our group-size halo (M200 = 1013 M) is quenched by a strong jet episode triggered by a cooling flow, and it is kept quenched by a low-power jet fed from hot halo accretion. In more massive systems (M200 ≳ 1014 M), hot halo accretion is insufficient to quench the galaxies, or to keep them quenched after the first cooling episode. These galaxies experience multiple episodes of gas cooling, star formation, and jet feedback. In the most massive galaxy cluster that we simulate (M200 = 1015 M), we find peak cold gas masses of 1010 M and peak star formation rates of a few times 100 $\mathrm{M}_\odot \,\, \mathrm{yr}^{-1}$. These values are achieved during strong cooling flows, which also trigger the strongest jets with peak powers of 1047$\mathrm{erg}\, \mathrm{s}^{-1}$. These jets subsequently shut off the cooling flows and any associated star formation. Jet-inflated bubbles draw out low-entropy gas that subsequently forms dense cooling filaments in their wakes, as seen in observations.

25) Immediate Origin of the Moon as a Post-impact Satellite

  Kegerreis, J et al., ApJL, 2022 , vol. 937 , issue 2   (citations: 4)


The Moon is traditionally thought to have coalesced from the debris ejected by a giant impact onto the early Earth. However, such models struggle to explain the similar isotopic compositions of Earth and lunar rocks at the same time as the system's angular momentum, and the details of potential impact scenarios are hotly debated. Above a high resolution threshold for simulations, we find that giant impacts can immediately place a satellite with similar mass and iron content to the Moon into orbit far outside Earth's Roche limit. Even satellites that initially pass within the Roche limit can reliably and predictably survive, by being partially stripped and then torqued onto wider, stable orbits. Furthermore, the outer layers of these directly formed satellites are molten over cooler interiors and are composed of around 60% proto-Earth material. This could alleviate the tension between the Moon's Earth-like isotopic composition and the different signature expected for the impactor. Immediate formation opens up new options for the Moon's early orbit and evolution, including the possibility of a highly tilted orbit to explain the lunar inclination, and offers a simpler, single-stage scenario for the origin of the Moon.

24) The importance of black hole repositioning for galaxy formation simulations

  Bahe, Y et al., MNRAS, 2022 , vol. 516 , issue 1   (citations: 20)


Active galactic nucleus (AGN) feedback from accreting supermassive black holes (SMBHs) is an essential ingredient of galaxy formation simulations. The orbital evolution of SMBHs is affected by dynamical friction that cannot be predicted self-consistently by contemporary simulations of galaxy formation in representative volumes. Instead, such simulations typically use a simple 'repositioning' of SMBHs, but the effects of this approach on SMBH and galaxy properties have not yet been investigated systematically. Based on a suite of smoothed particle hydrodynamics simulations with the SWIFT code and a Bondi-Hoyle-Lyttleton sub-grid gas accretion model, we investigate the impact of repositioning on SMBH growth and on other baryonic components through AGN feedback. Across at least a factor ~1000 in mass resolution, SMBH repositioning (or an equivalent approach) is a necessary prerequisite for AGN feedback; without it, black hole growth is negligible. Limiting the effective repositioning speed to ≲10 km s-1 delays the onset of AGN feedback and severely limits its impact on stellar mass growth in the centre of massive galaxies. Repositioning has three direct physical consequences. It promotes SMBH mergers and thus accelerates their initial growth. In addition, it raises the peak density of the ambient gas and reduces the SMBH velocity relative to it, giving a combined boost to the accretion rate that can reach many orders of magnitude. Our results suggest that a more sophisticated and/or better calibrated treatment of SMBH repositioning is a critical step towards more predictive galaxy formation simulations.

23) The interplay between AGN feedback and precipitation of the intracluster medium in simulations of galaxy groups and clusters

  Nobels, F et al., MNRAS, 2022 , vol. 515 , issue 4   (citations: 11)


Using high-resolution hydrodynamical simulations of idealized galaxy clusters, we study the interaction between the brightest cluster galaxy, its supermassive black hole (BH), and the intracluster medium (ICM). We create initial conditions for which the ICM is in hydrostatic equilibrium within the gravitational potential from the galaxy and an NFW dark matter halo. Two free parameters associated with the thermodynamic profiles determine the cluster gas fraction and the central temperature, where the latter can be used to create cool-core or non-cool-core systems. Our simulations include radiative cooling, star formation, BH accretion, and stellar and active galactic nucleus (AGN) feedback. Even though the energy of AGN feedback is injected thermally and isotropically, it leads to anisotropic outflows and buoyantly rising bubbles. We find that the BH accretion rate (BHAR) is highly variable and only correlates strongly with the star formation rate (SFR) and the ICM when it is averaged over more than $1~\rm Myr$. We generally find good agreement with the theoretical precipitation framework. In $10^{13}~\rm M_\odot$ haloes, AGN feedback quenches the central galaxy and converts cool-core systems into non-cool-core systems. In contrast, higher mass, cool-core clusters evolve cyclically. Episodes of high BHAR raise the entropy of the ICM out to the radius, where the ratio of the cooling time and the local dynamical time tcool/tdyn > 10, thus suppressing condensation and, after a delay, the BHAR. The corresponding reduction in AGN feedback allows the ICM to cool and become unstable to precipitation, thus initiating a new episode of high SFR and BHAR.

22) Continuous Simulation Data Stream: A dynamical timescale-dependent output scheme for simulations

  Hausammann, L et al., A&C, 2022 , vol. 41   (citations: 1)


Exa-scale simulations are on the horizon but almost no new design for the output has been proposed in recent years. In simulations using individual time steps, the traditional snapshots are over resolving particles/cells with large time steps and are under resolving the particles/cells with short time steps. Therefore, they are unable to follow fast events and use efficiently the storage space. The Continuous Simulation Data Stream (CSDS) is designed to decrease this space while providing an accurate state of the simulation at any time. It takes advantage of the individual time step to ensure the same relative accuracy for all the particles. The outputs consist of a single file representing the full evolution of the simulation. Within this file, the particles are written independently and at their own frequency. Through the interpolation of the records, the state of the simulation can be recovered at any point in time. In this paper, we show that the CSDS can reduce the storage space by 2.76x for the same accuracy than snapshots or increase the accuracy by 67.8x for the same storage space whilst retaining an acceptable reading speed for analysis. By using interpolation between records, the CSDS provides the state of the simulation, with a high accuracy, at any time. This should largely improve the analysis of fast events such as supernovae and simplify the construction of light-cone outputs.

21) The DESI N-body simulation project - I. Testing the robustness of simulations for the DESI dark time survey

  Grove, C et al., MNRAS, 2022 , vol. 515 , issue 2   (citations: 16)


Analysis of large galaxy surveys requires confidence in the robustness of numerical simulation methods. The simulations are used to construct mock galaxy catalogues to validate data analysis pipelines and identify potential systematics. We compare three N-body simulation codes, ABACUS, GADGET-2, and SWIFT, to investigate the regimes in which their results agree. We run N-body simulations at three different mass resolutions, 6.25 × 108, 2.11 × 109, and 5.00 × 109 h-1 M, matching phases to reduce the noise within the comparisons. We find systematic errors in the halo clustering between different codes are smaller than the Dark Energy Spectroscopic Instrument (DESI) statistical error for $s\ \gt\ 20\ h^{-1}$ Mpc in the correlation function in redshift space. Through the resolution comparison we find that simulations run with a mass resolution of 2.1 × 109 h-1 M are sufficiently converged for systematic effects in the halo clustering to be smaller than the DESI statistical error at scales larger than $20\ h^{-1}$ Mpc. These findings show that the simulations are robust for extracting cosmological information from large scales which is the key goal of the DESI survey. Comparing matter power spectra, we find the codes agree to within 1 per cent for k ≤ 10 h Mpc-1. We also run a comparison of three initial condition generation codes and find good agreement. In addition, we include a quasi-N-body code, FastPM, since we plan use it for certain DESI analyses. The impact of the halo definition and galaxy-halo relation will be presented in a follow-up study.

20) The importance of the way in which supernova energy is distributed around young stellar populations in simulations of galaxies

  Chaikin, E et al., MNRAS, 2022 , vol. 514 , issue 1   (citations: 12)


Supernova (SN) feedback plays a crucial role in simulations of galaxy formation. Because blast waves from individual SNe occur on scales that remain unresolved in modern cosmological simulations, SN feedback must be implemented as a subgrid model. Differences in the manner in which SN energy is coupled to the local interstellar medium and in which excessive radiative losses are prevented have resulted in a zoo of models used by different groups. However, the importance of the selection of resolution elements around young stellar particles for SN feedback has largely been overlooked. In this work, we examine various selection methods using the smoothed particle hydrodynamics code SWIFT. We run a suite of isolated disc galaxy simulations of a Milky Way-mass galaxy and small cosmological volumes, all with the thermal stochastic SN feedback model used in the EAGLE simulations. We complement the original mass-weighted neighbour selection with a novel algorithm guaranteeing that the SN energy distribution is as close to isotropic as possible. Additionally, we consider algorithms where the energy is injected into the closest, least dense, or most dense neighbour. We show that different neighbour-selection strategies cause significant variations in star formation rates, gas densities, wind mass-loading factors, and galaxy morphology. The isotropic method results in more efficient feedback than the conventional mass-weighted selection. We conclude that the manner in which the feedback energy is distributed among the resolution elements surrounding a feedback event is as important as changing the amount of energy by factors of a few.

19) SIBELIUS-DARK: a galaxy catalogue of the local volume from a constrained realization simulation

  McAlpine, S et al., MNRAS, 2022 , vol. 512 , issue 4   (citations: 17)


We present SIBELIUS-DARK, a constrained realization simulation of the local volume to a distance of 200 Mpc from the Milky Way. SIBELIUS-DARK is the first study of the 'Simulations Beyond The Local Universe' (SIBELIUS) project, which has the goal of embedding a model Local Group-like system within the correct cosmic environment. The simulation is dark-matter-only, with the galaxy population calculated using the semi-analytic model of galaxy formation, GALFORM. We demonstrate that the large-scale structure that emerges from the SIBELIUS constrained initial conditions matches well the observational data. The inferred galaxy population of SIBELIUS-DARK also match well the observational data, both statistically for the whole volume and on an object-by-object basis for the most massive clusters. For example, the K-band number counts across the whole sky, and when divided between the northern and southern Galactic hemispheres, are well reproduced by SIBELIUS-DARK. We find that the local volume is somewhat unusual in the wider context of ΛCDM: it contains an abnormally high number of supermassive clusters, as well as an overall large-scale underdensity at the level of ≈5 per cent relative to the cosmic mean. However, whilst rare, the extent of these peculiarities does not significantly challenge the ΛCDM model. SIBELIUS-DARK is the most comprehensive constrained realization simulation of the local volume to date, and with this paper we publicly release the halo and galaxy catalogues at z = 0, which we hope will be useful to the wider astronomy community.

18) Simulations of 60Fe entrained in ejecta from a near-Earth supernova: effects of observer motion

  Chaikin, E et al., MNRAS, 2022 , vol. 512 , issue 1   (citations: 6)


Recent studies have shown that live (not decayed) radioactive 60Fe is present in deep-ocean samples, Antarctic snow, lunar regolith, and cosmic rays. 60Fe represents supernova (SN) ejecta deposited in the Solar system around $3 \, \rm Myr$ ago, and recently an earlier pulse ${\approx}7 \ \rm Myr$ ago has been found. These data point to one or multiple near-Earth SN explosions that presumably participated in the formation of the Local Bubble. We explore this theory using 3D high-resolution smooth-particle hydrodynamical simulations of isolated SNe with ejecta tracers in a uniform interstellar medium (ISM). The simulation allows us to trace the SN ejecta in gas form and those eject in dust grains that are entrained with the gas. We consider two cases of diffused ejecta: when the ejecta are well-mixed in the shock and when they are not. In the latter case, we find that these ejecta remain far behind the forward shock, limiting the distance to which entrained ejecta can be delivered to ≈100 pc in an ISM with $n_\mathrm{H}=0.1\,\, \rm cm^{-3}$ mean hydrogen density. We show that the intensity and the duration of 60Fe accretion depend on the ISM density and the trajectory of the Solar system. Furthermore, we show the possibility of reproducing the two observed peaks in 60Fe concentration with this model by assuming two linear trajectories for the Solar system with 30-km s-1 velocity. The fact that we can reproduce the two observed peaks further supports the theory that the 60Fe signal was originated from near-Earth SNe.

17) Fast full N-body simulations of generic modified gravity: conformal coupling models

  Ruan, C et al., JCAP, 2022 , vol. 2022 , issue 5   (citations: 15)


We present MG-GLAM, a code developed for the very fast production of full N-body cosmological simulations in modified gravity (MG) models. We describe the implementation, numerical tests and first results of a large suite of cosmological simulations for three classes of MG models with conformal coupling terms: the f(R) gravity, symmetron and coupled quintessence models. Derived from the parallel particle-mesh code GLAM, MG-GLAM incorporates an efficient multigrid relaxation technique to solve the characteristic nonlinear partial differential equations of these models. For f(R) gravity, we have included new variants to diversify the model behaviour, and we have tailored the relaxation algorithms to these to maintain high computational efficiency. In a companion paper, we describe versions of this code developed for derivative coupling MG models, including the Vainshtein- and K-mouflage-type models. MG-GLAM can model the prototypes for most MG models of interest, and is broad and versatile. The code is highly optimised, with a tremendous speedup of a factor of more than a hundred compared with earlier N-body codes, while still giving accurate predictions of the matter power spectrum and dark matter halo abundance. MG-GLAM is ideal for the generation of large numbers of MG simulations that can be used in the construction of mock galaxy catalogues and the production of accurate emulators for ongoing and future galaxy surveys.

16) SPHENIX: smoothed particle hydrodynamics for the next generation of galaxy formation simulations

  Borrow, J et al., MNRAS, 2022 , vol. 511 , issue 2   (citations: 27)


Smoothed particle hydrodynamics (SPH) is a ubiquitous numerical method for solving the fluid equations, and is prized for its conservation properties, natural adaptivity, and simplicity. We introduce the SPHENIX SPH scheme, which was designed with three key goals in mind: to work well with sub-grid physics modules that inject energy, be highly computationally efficient (both in terms of compute and memory), and to be Lagrangian. SPHENIX uses a Density-Energy equation of motion, along with a variable artificial viscosity and conduction, including limiters designed to work with common sub-grid models of galaxy formation. In particular, we present and test a novel limiter that prevents conduction across shocks, preventing spurious radiative losses in feedback events. SPHENIX is shown to solve many difficult test problems for traditional SPH, including fluid mixing and vorticity conservation, and it is shown to produce convergent behaviour in all tests where this is appropriate. Crucially, we use the same parameters within SPHENIX for the various switches throughout, to demonstrate the performance of the scheme as it would be used in production simulations. SPHENIX is the new default scheme in the SWIFT cosmological simulation code and is available open source.

15) swift-emulator: A Python package for emulation of simulated scaling relations

  Kugel, R et al., JOSS, 2022 , vol. 7 , issue 72   (citations: 6)

    No abstract

14) Snowmass2021 Cosmic Frontier White Paper: High Density Galaxy Clustering in the Regime of Cosmic Acceleration

  Dawson, K et al., arXiv, 2022   (citations: 11)


Joint studies of imaging and spectroscopic samples, informed by theory and simulations, offer the potential for comprehensive tests of the cosmological model over redshifts z<1.5. Spectroscopic galaxy samples at these redshifts can be increased beyond the planned Dark Energy Spectroscopic Instrument (DESI) program by at least an order of magnitude, thus offering significantly more constraining power for these joint studies. Spectroscopic observations of these galaxies in the latter half of the 2020's and beyond would leverage the theory and simulation effort in this regime. In turn, these high density observations will allow enhanced tests of dark energy, physics beyond the standard model, and neutrino masses that will greatly exceed what is currently possible. Here, we present a coordinated program of simulations, theoretical modeling, and future spectroscopy that would enable precise cosmological studies in the accelerating epoch where the effects of dark energy are most apparent.

13) Fast full N-body simulations of generic modified gravity: derivative coupling models

  Hernandez-Aguayo, C et al., JCAP, 2022 , vol. 2022 , issue 1   (citations: 13)


We present MG-GLAM, a code developed for the very fast production of full N-body cosmological simulations in modified gravity (MG) models. We describe the implementation, numerical tests and first results of a large suite of cosmological simulations for two broad classes of MG models with derivative coupling terms - the Vainshtein- and Kmouflage-type models - which respectively features the Vainshtein and Kmouflage screening mechanism. Derived from the parallel particle-mesh code GLAM, MG-GLAM incorporates an efficient multigrid relaxation technique to solve the characteristic nonlinear partial differential equations of these models. For Kmouflage, we have proposed a new algorithm for the relaxation solver, and run the first simulations of the model to understand its cosmological behaviour. In a companion paper, we describe versions of this code developed for conformally-coupled MG models, including several variants of f(R) gravity, the symmetron model and coupled quintessence. Altogether, MG-GLAM has so far implemented the prototypes for most MG models of interest, and is broad and versatile. The code is highly optimised, with a tremendous (over two orders of magnitude) speedup when comparing its running time with earlier N-body codes, while still giving accurate predictions of the matter power spectrum and dark matter halo abundance. MG-GLAM is ideal for the generation of large numbers of MG simulations that can be used in the construction of mock galaxy catalogues and accurate emulators for ongoing and future galaxy surveys.

12) Smoothed particle radiation hydrodynamics: two-moment method with local Eddington tensor closure

  Chan, T et al., MNRAS, 2021 , vol. 505 , issue 4   (citations: 10)


We present a new smoothed particle hydrodynamics-radiative transfer method (SPH-M1RT) that is coupled dynamically with SPH. We implement it in the (task-based parallel) SWIFT galaxy simulation code but it can be straightforwardly implemented in other SPH codes. Our moment-based method simultaneously solves the radiation energy and flux equations in SPH, making it adaptive in space and time. We modify the M1 closure relation to stabilize radiation fronts in the optically thin limit. We also introduce anisotropic artificial viscosity and high-order artificial diffusion schemes, which allow the code to handle radiation transport accurately in both the optically thin and optically thick regimes. Non-equilibrium thermochemistry is solved using a semi-implicit sub-cycling technique. The computational cost of our method is independent of the number of sources and can be lowered further by using the reduced speed-of-light approximation. We demonstrate the robustness of our method by applying it to a set of standard tests from the cosmological radiative transfer comparison project of Iliev et al. The SPH-M1RT scheme is well-suited for modelling situations in which numerous sources emit ionizing radiation, such as cosmological simulations of galaxy formation or simulations of the interstellar medium.

11) Nyx: A Massively Parallel AMR Code for Computational Cosmology

  Sexton, J et al., JOSS, 2021 , vol. 6 , issue 63   (citations: 4)

    No abstract

10) Higher order initial conditions for mixed baryon-CDM simulations

  Hahn, O et al., MNRAS, 2021 , vol. 503 , issue 1   (citations: 25)


We present a novel approach to generate higher order initial conditions (ICs) for cosmological simulations that take into account the distinct evolution of baryons and dark matter. We focus on the numerical implementation and the validation of its performance, based on both collisionless N-body simulations and full hydrodynamic Eulerian and Lagrangian simulations. We improve in various ways over previous approaches that were limited to first-order Lagrangian perturbation theory (LPT). Specifically, we (1) generalize nth-order LPT to multifluid systems, allowing 2LPT or 3LPT ICs for two-fluid simulations, (2) employ a novel propagator perturbation theory to set up ICs for Eulerian codes that are fully consistent with 1LPT or 2LPT, (3) demonstrate that our ICs resolve previous problems of two-fluid simulations by using variations in particle masses that eliminate spurious deviations from expected perturbative results, (4) show that the improvements achieved by going to higher order PT are comparable to those seen for single-fluid ICs, and (5) demonstrate the excellent (i.e. few per cent level) agreement between Eulerian and Lagrangian simulations, once high-quality initial conditions are used. The rigorous development of the underlying perturbation theory is presented in a companion paper. All presented algorithms are implemented in the MONOFONIC MUSIC-2 package that we make publicly available.

9) The imprint of dark subhaloes on the circumgalactic medium

  McCarthy, I et al., MNRAS, 2020 , vol. 499 , issue 3   (citations: 1)


The standard model of cosmology, the Λ cold dark matter (ΛCDM) model, robustly predicts the existence of a multitude of dark matter 'subhaloes' around galaxies like the Milky Way. A wide variety of observations have been proposed to look for the gravitational effects such subhaloes would induce in observable matter. Most of these approaches pertain to the stellar or cool gaseous phases of matter. Here we propose a new approach, which is to search for the perturbations that such dark subhaloes would source in the warm/hot circumgalactic medium (CGM) around normal galaxies. With a combination of analytic theory, carefully controlled high-resolution idealized simulations, and full cosmological hydrodynamical simulations (the ARTEMIS simulations), we calculate the expected signal and how it depends on important physical parameters (subhalo mass, CGM temperature, and relative velocity). We find that dark subhaloes enhance both the local CGM temperature and density and, therefore, also the pressure. For the pressure and density, the fluctuations can vary in magnitude from tens of per cent (for subhaloes with Msub = 1010 M) to a few per cent (for subhaloes with Msub = 108 M), although this depends strongly on the CGM temperature. The subhaloes also induce fluctuations in the velocity field ranging in magnitude from a few km s-1 up to 25 km s-1. We propose that X-ray, Sunyaev-Zel'dovich effect, radio dispersion measure, and quasar absorption line observations can be used to measure these fluctuations and place constraints on the abundance and distribution of dark subhaloes, thereby placing constraints on the nature of dark matter.

8) A versatile smoothed particle hydrodynamics code for graphic cards

  Schafer, C et al., A&C, 2020 , vol. 33   (citations: 7)


We present the second release of the now open source smoothed particle hydrodynamics code miluphcuda. The code is designed to run on Nvidia CUDA capable devices. It handles one to three dimensional problems and includes modules to solve the equations for viscid and inviscid hydrodynamical flows, the equations of continuum mechanics using SPH, and self-gravity with a Barnes-Hut tree. The covered material models include different porosity and plasticity models. Several equations of states, especially for impact physics, are implemented. The basic ideas of the numerical scheme are presented, the usage of the code is explained and its versatility is shown by means of different applications. The code is hereby publicly available.

7) CMACIONIZE 2.0: a novel task-based approach to Monte Carlo radiation transfer

  Vandenbroucke, B et al., A&A, 2020 , vol. 641   (citations: 2)


Context. Monte Carlo radiative transfer (MCRT) is a widely used technique to model the interaction between radiation and a medium. It plays an important role in astrophysical modelling and when these models are compared with observations.
Aims: We present a novel approach to MCRT that addresses the challenging memory-access patterns of traditional MCRT algorithms, which prevent an optimal performance of MCRT simulations on modern hardware with a complex memory architecture.
Methods: We reformulated the MCRT photon-packet life cycle as a task-based algorithm, whereby the computation is broken down into small tasks that are executed concurrently. Photon packets are stored in intermediate buffers, and tasks propagate photon packets through small parts of the computational domain, moving them from one buffer to another in the process.
Results: Using the implementation of the new algorithm in the photoionization MCRT code CMACIONIZE 2.0, we show that the decomposition of the MCRT grid into small parts leads to a significant performance gain during the photon-packet propagation phase, which constitutes the bulk of an MCRT algorithm because memory caches are used more efficiently. Our new algorithm is faster by a factor 2 to 4 than an equivalent traditional algorithm and shows good strong scaling up to 30 threads. We briefly discuss adjustments to our new algorithm and extensions to other astrophysical MCRT applications.
Conclusions: We show that optimising the memory access patterns of a memory-bound algorithm such as MCRT can yield significant performance gains. The source code of CMACIONIZE 2.0 is hosted at https://github.com/bwvdnbro/CMacIonize

6) swiftsimio: A Python library for reading SWIFT data

  Borrow, J et al., JOSS, 2020 , vol. 5 , issue 52   (citations: 32)

    No abstract

5) Atmospheric Erosion by Giant Impacts onto Terrestrial Planets

  Kegerreis, J et al., ApJ, 2020 , vol. 897 , issue 2   (citations: 21)


We examine the mechanisms by which the atmosphere can be eroded by giant impacts onto Earth-like planets with thin atmospheres, using 3D smoothed particle hydrodynamics simulations with sufficient resolution to directly model the fate of low-mass atmospheres. We present a simple scaling law to estimate the fraction lost for any impact angle and speed in this regime. In the canonical Moon-forming impact, only around 10% of the atmosphere would have been lost from the immediate effects of the collision. There is a gradual transition from removing almost none to almost all of the atmosphere for a grazing impact as it becomes more head-on or increases in speed, including complex, nonmonotonic behavior at low impact angles. In contrast, for head-on impacts, a slightly greater speed can suddenly remove much more atmosphere. Our results broadly agree with the application of 1D models of local atmosphere loss to the ground speeds measured directly from our simulations. However, previous analytical models of shock-wave propagation from an idealized point-mass impact significantly underestimate the ground speeds and hence the total erosion. The strong dependence on impact angle and the interplay of multiple nonlinear and asymmetrical loss mechanisms highlight the need for 3D simulations in order to make realistic predictions.

4) Honing and proofing Astrophysical codes on the road to Exascale. Experiences from code modernization on many-core systems

  Cielo, S et al., arXiv, 2020   (citations: 0)


The complexity of modern and upcoming computing architectures poses severe challenges for code developers and application specialists, and forces them to expose the highest possible degree of parallelism, in order to make the best use of the available hardware. The Intel$^{(R)}$ Xeon Phi$^{(TM)}$ of second generation (code-named Knights Landing, henceforth KNL) is the latest many-core system, which implements several interesting hardware features like for example a large number of cores per node (up to 72), the 512 bits-wide vector registers and the high-bandwidth memory. The unique features of KNL make this platform a powerful testbed for modern HPC applications. The performance of codes on KNL is therefore a useful proxy of their readiness for future architectures. In this work we describe the lessons learnt during the optimisation of the widely used codes for computational astrophysics P-Gadget-3, Flash and Echo. Moreover, we present results for the visualisation and analysis tools VisIt and yt. These examples show that modern architectures benefit from code optimisation at different levels, even more than traditional multi-core systems. However, the level of modernisation of typical community codes still needs improvements, for them to fully utilise resources of novel architectures.

3) Hunting for galaxies and halos in simulations with VELOCIraptor

  Elahi, P et al., PASA, 2019 , vol. 36   (citations: 68)


We present VELOCIraptor, a massively parallel galaxy/(sub)halo finder that is also capable of robustly identifying tidally disrupted objects and separate stellar halos from galaxies. The code is written in C++11, use the Message Passing Interface (MPI) and OpenMP Application Programming Interface (API) for parallelisation, and includes python tools to read/manipulate the data products produced. We demonstrate the power of the VELOCIraptor (sub)halo finder, showing how it can identify subhalos deep within the host that have negligible density contrasts to their parent halo. We find a subhalo mass-radial distance dependence: large subhalos with mass ratios of ≳10-2 are more common in the central regions than smaller subhalos, a result of dynamical friction and low tidal mass loss rates. This dependence is completely absent in (sub)halo finders in common use, which generally search for substructure in configuration space, yet is present in codes that track particles belonging to halos as they fall into other halos, such as hbt+. VELOCIraptor largely reproduces the dependence seen without tracking, finding a similar radial dependence to hbt+ in well-resolved halos from our limited resolution fiducial simulation.

2) A momentum conserving $N$-body scheme with individual timesteps

  Zhu, Q, arXiv, 2017   (citations: 2)


$N$-body simulations study the dynamics of $N$ particles under the influence of mutual long-distant forces such as gravity. In practice, $N$-body codes will violate Newton's third law if they use either an approximate Poisson solver or individual timesteps. In this study, we construct a novel $N$-body scheme by combining a fast multipole method (FMM) based Poisson solver and a time integrator using a hierarchical Hamiltonian splitting (HHS) technique. We test our implementation for collision-less systems using several problems in galactic dynamics. As a result of the momentum conserving nature of these two key components, the new $N$-body scheme is also momentum conserving. Moreover, we can fully utilize the $\mathcal O(\textit N)$ complexity of FMM with the integrator. With the restored force symmetry, we can improve both angular momentum conservation and energy conservation substantially. The new scheme will be suitable for many applications in galactic dynamics and structure formation. Our implementation, in the code Taichi, is publicly available at https://bitbucket.org/qirong_zhu/taichi_public/.

1) SpECTRE: A task-based discontinuous Galerkin code for relativistic astrophysics

  Kidder, L et al., JCoPh, 2017 , vol. 335   (citations: 73)


We introduce a new relativistic astrophysics code, SpECTRE, that combines a discontinuous Galerkin method with a task-based parallelism model. SpECTRE's goal is to achieve more accurate solutions for challenging relativistic astrophysics problems such as core-collapse supernovae and binary neutron star mergers. The robustness of the discontinuous Galerkin method allows for the use of high-resolution shock capturing methods in regions where (relativistic) shocks are found, while exploiting high-order accuracy in smooth regions. A task-based parallelism model allows efficient use of the largest supercomputers for problems with a heterogeneous workload over disparate spatial and temporal scales. We argue that the locality and algorithmic structure of discontinuous Galerkin methods will exhibit good scalability within a task-based parallelism framework. We demonstrate the code on a wide variety of challenging benchmark problems in (non)-relativistic (magneto)-hydrodynamics. We demonstrate the code's scalability including its strong scaling on the NCSA Blue Waters supercomputer up to the machine's full capacity of 22 , 380 nodes using 671 , 400 threads.