Direct N-body Simulations Research Articles

Origin of high dark remnant fractions in Milky Way Globular Clusters: The crucial role of initial black hole retention

Abstract Comparing the dynamical and stellar masses of Milky Way (MW) globular clusters (GCs) reveals a discrepancy exceeding a factor of two. Since this substantial invisible mass is concentrated in the cluster centre, it is attributed to stellar remnants. The majority of mass in remnants consists of white dwarfs (WDs). Allocating over half of a GC’s current mass to WDs could significantly restrict the dynamical evolution scenarios governing stellar clusters. As the most massive stars in GCs, black holes (BHs) exert a substantial effect on the escape rate of lower mass stars, such as WDs. This paper aims to identify which scenarios of BH natal kicks can accurately reproduce the notable dark remnant fraction observed in MW GCs. We compare the observed remnant fraction of MW GCs with a comprehensive grid of direct N-body simulations while adjusting the natal kick received by BHs. Our results reveal that simulations employing low natal kicks to BHs are the only ones capable of mirroring the remnant fraction of MW GCs. According to the Spitzer instability, the presence of a BH population prompts the formation of a BH sub-system (BHSub) at the centre of a star cluster. The BHSub serves as an energetic power plant, continually releasing kinetic energy through few-body encounters between single and binary BHs, and transferring the generated energy to the entire stellar population. This energy induces a significant difference in the ejection rate of stellar remnants and luminous stars, ultimately increasing the fraction of dark remnants within the star cluster.

The Potential for Long-lived Intermediate-mass Black Hole Binaries in the Lowest Density Dwarf Galaxies

Intermediate-mass black hole (IMBH) mergers with masses 104–106 M ⊙ are expected to produce gravitational waves detectable by the Laser Interferometer Space Antenna (LISA) with high signal-to-noise ratios from the present day to cosmic dawn. IMBH mergers are expected to take place within dwarf galaxies; however, the dynamics, timescales, and effect on their hosts are largely unexplored. In a previous study, we examined how IMBHs would pair and merge within nucleated dwarf galaxies. IMBHs in nucleated hosts evolve very efficiently, forming a binary system and coalescing within a few hundred million years. Although the fraction of dwarf galaxies (107 M ⊙ ≤ M ⋆ ≤ 1010 M ⊙) hosting nuclear star clusters is between 60% and 100%, this fraction drops to 20%–70% for lower-mass dwarfs (M ⋆ ≈ 107 M ⊙), with the largest drop in low-density environments. Here, we extend our previous study by performing direct N-body simulations to explore the dynamics and evolution of IMBHs within nonnucleated dwarf galaxies, under the assumption that IMBHs exist within these dwarfs. To our surprise, none of the IMBHs in our simulation suite merge within a Hubble time, despite many attaining high eccentricities e ∼ 0.7–0.95. We conclude that extremely low stellar density environments in the centers of nonnucleated dwarfs do not provide an ample supply of stars to interact with an IMBH binary, resulting in its stalling, in spite of triaxiality and high eccentricity, common means to drive a binary to coalescence. Our findings underline the importance of considering all detailed host properties to predict IMBH merger rates for LISA.

A Systematic Analysis of Star Cluster Disruption by Tidal Shocks. II. Predicting Star Cluster Dissolution Rates from a Time-series Analysis of Their Tidal Histories

Most of the dynamical mass loss from star clusters is thought to be caused by the time variability of the tidal field (“tidal shocks”). Systematic studies of tidal shocks have been hampered by the fact that each tidal history is unique, implying both a reproducibility and a generalization problem. Here we address these issues by investigating how star cluster evolution depends on the statistical properties of its tidal history. We run a large suite of direct N-body simulations of clusters with tidal histories generated from power spectra of a given slope and with different normalizations, which determine the timescales and amplitudes of the shocks, respectively. At fixed normalization (i.e., the same median tidal field strength), the dissolution timescale is nearly independent of the power spectrum slope. However, the dispersion in dissolution timescales, obtained by repeating simulations for different realizations of statistically identical tidal histories, increases with the power spectrum slope. This result means that clusters experiencing high-frequency shocks have more similar mass-loss histories than clusters experiencing low-frequency shocks. The density–mass relationship of the simulated clusters follows a power law with slope between 1.08 and 1.45, except for the lowest normalizations (for which clusters effectively evolve in a static tidal field). Our findings suggest that star cluster evolution can be described statistically from a time-series analysis of its tidal history, which is an important simplification for describing the evolution of the star cluster population during galaxy formation and evolution.

Evolution of the mass-radius relation of expanding very young star clusters

The initial mass–radius relation of embedded star clusters is an essential boundary condition for understanding the evolution of embedded clusters in which stars form to their release into the galactic field via an open star cluster phase. The initial mass–radius relation of embedded clusters deduced by Marks & Kroupa (2012, A&A, 543, A8) is significantly different from the relation suggested by Pfalzner et al. (2016, A&A, 586, A68). Here, we use direct N-body simulations to model the early expansion of embedded clusters after the expulsion of their residual gas. The observationally deduced radii of clusters up to a few million years old, compiled from various sources, are well fitted by N-body models, implying that these observed very young clusters are most likely in an expanding state. We show that the mass–radius relation of Pfalzner et al. (2016) reflects the expansion of embedded clusters following the initial mass–radius relation of Marks & Kroupa (2012). We also suggest that even the embedded clusters in ATLASGAL clumps with HII regions are probably already in expansion. All the clusters collected here from different observations show a mass-radius relation with a similar slope, which may indicate that all clusters were born with a profile resembling that of the Plummer phase-space distribution function.

On the Dynamical Evolution of the Asteroid Belt in a Massive Star–Neutron Star Binary

Some fast radio bursts (FRBs) exhibit repetitive behaviors, and their origins remain enigmatic. It has been argued that repeating FRBs could be produced by the interaction between a neutron star and an asteroid belt. Here, we consider the systems in which an asteroid belt dwells around a massive star, while a neutron star, as a companion of the massive star, interacts with the belt through gravitational force. Various orbital configurations are assumed for the system. Direct N-body simulations are performed to investigate the dynamical evolution of the asteroids' belt. It is found that a larger orbital eccentricity of the neutron star will destroy the belt more quickly, with a large number of asteroids being scattered out of the system. A low inclination not only suppresses the collisions but also inhibits the ejection rate at early stages. However, highly inclined systems may undergo strong oscillations, resulting in the Kozai–Lidov instabilities. Among the various configurations, a clear periodicity is observed in the collision events for the case with an orbital eccentricity of 0.7 and mutual inclination of 0◦. It is found that such a periodicity can be sustained for at least eight neutron star orbital periods, supporting this mechanism as a possible explanation for periodically repeating FRBs. Our studies also suggest that the active stage of these kinds of FRB sources should be limited, since the asteroid belt would finally be destroyed by the neutron star after multiple passages.

Dynamical and Secular Stability of Mutually Inclined Planetary Systems

Multiple analytical, semi-analytical, and empirical stability criteria have been derived in the literature for two-planet systems. But, the dependence of the stability limit on the initial mutual inclination between the inner and outer orbits is not well modeled by previous stability criteria. Here, we derive a semi-analytical stability criteria for two-planet systems, at arbitrary inclinations, in which the inner planet is a test particle. Using perturbation theory we calculate the characteristic fractional change in the semimajor axis of the inner binary β = δ a 1/a 1 caused by perturbations from the companion. Stability criteria can be derived by setting a threshold on β. Focusing initially on circular orbits, we derive an analytical expression for β for coplanar prograde and retrograde orbits. For noncoplanar configurations, we evaluate a semi-analytical expression. We then generalize to orbits with arbitrary eccentricities and account for the secular effects. Our analytical and semi-analytical results are in excellent agreement with direct N-body simulations. In addition, we show that contours of β ∼ 0.01 can serve as criteria for stability. More specifically, we show that (1) retrograde orbits are generally more stable than prograde ones; (2) systems with intermediate mutual inclination are less stable due to von Ziepel–Lidov–Kozai (vZLK) dynamics; and (3) mean motion resonances (MMRs) can stabilize intermediate inclination secularly unstable regions in phase space, by quenching vZLK secular processes (4) MMRs can destabilize some of the dynamically stable regions. We also point out that these stability criteria can be used to constrain the orbital properties of observed systems and their age.

Extensive Pollution of Uranus and Neptune’s Atmospheres by Upsweep of Icy Material during the Nice Model Migration

In the Nice model of Solar System formation, Uranus and Neptune undergo an orbital upheaval, sweeping through a planetesimal disk. The region of the disk from which material is accreted by the ice giants during this phase of their evolution has not previously been identified. We perform direct N-body orbital simulations of the four giant planets to determine the amount and origin of solid accretion during this orbital upheaval. We find that the ice giants undergo an extreme bombardment event, with collision rates as high as ∼3 per hour assuming km-sized planetesimals, increasing the total planet mass by up to ∼0.35%. In all cases, the initially outermost ice giant experiences the largest total enhancement. We determine that, for some plausible planetesimal properties, the resulting atmospheric enrichment could potentially produce sufficient latent heat to alter the planetary cooling timescale according to existing models. Our findings suggest that substantial accretion during this phase of planetary evolution may have been sufficient to impact the atmospheric composition and thermal evolution of the ice giants, motivating future work on the fate of deposited solid material.

The impact of a top-heavy IMF on the formation and evolution of dark star clusters

ABSTRACT The Spitzer instability leads to the formation of a black hole subsystem (BHSub) at the centre of a star cluster providing energy to luminous stars (LSs) and increasing their rate of evaporation. When the self-depletion time of the BHSub exceeds the evaporation time of the LSs, a dark star cluster (DSC) will appear. Using the nbody7 code, we performed a comprehensive set of direct N-body simulations over a wide range of initial conditions to study the pure effect of the top-heaviness of the initial mass function (IMF) on the formation of the DSC phase. In the Galactic tidal field, top-heavy IMFs lead to the fast evaporation of LSs and the formation of DSCs. Therefore, DSCs can be present even in the outer region of the Milky Way (MW). To successfully transition to the DSC phase, the MW Globular Clusters (GCs) must possess an initial black hole (BH) mass fraction of $\widetilde {\mathit {M}}_\mathrm{BH}(0)\gt 0.05$. For star clusters with $\widetilde {\mathit {M}}_\mathrm{BH}(0)\gt 0.08$, the DSC phase will be created for any given initial density of the cluster and Galactocentric distance. The duration of the cluster’s lifetime spent in the DSC phase shows a negative (positive) correlation with the initial density, and Galactocentric distance of the star cluster if $\widetilde {\mathit {M}}_\mathrm{BH}(0)\le 0.12$ ($\widetilde {\mathit {M}}_\mathrm{BH}(0)\ge 0.15$). Considering the canonical IMF, it is unlikely for any MW GCs to enter the DSC phase. We discuss the BH retention fraction in view of the observed properties of the GCs of the MW.

Young Star Clusters Dominate the Production of Detached Black Hole–Star Binaries

The recent discovery of two detached black hole–star (BH–star) binaries from Gaia’s third data release has sparked interest in understanding the formation mechanisms of these systems. We investigate the formation of these systems by dynamical processes in young star clusters (SCs) and via isolated binary (IB) evolution, using a combination of direct N-body and population synthesis simulations. We find that dynamical formation in SCs is nearly 50 times more efficient per unit of star formation at producing BH–star binaries than IB evolution. We expand this analysis to the full Milky Way (MW) using a FIRE-2 hydrodynamical simulation of an MW-mass galaxy. Even assuming that only 10% of star formation goes into SCs, we find that approximately four out of every five BH–star systems are formed dynamically, and that the MW contains a total of ∼2 × 105 BH–star systems. Many of these dynamically formed systems have longer orbital periods, greater eccentricities, and greater black hole masses than their isolated counterparts. For binaries older than 100 Myr, we show that any detectable system with e ≳ 0.5 or M BH ≳ 10 M ⊙ can only be formed through dynamical processes. Our MW model predicts between 64 and 215 such detections from the complete DR4 Gaia catalog, with the majority of systems being dynamically formed in massive and metal-rich SCs. Finally, we compare our populations to the recently discovered Gaia BH1 and Gaia BH2, and conclude that the dynamical scenario is the most favorable formation pathway for both systems.

Dynamics of 2023 FW14, the second L4 Mars trojan, and a physical characterization using the 10.4 m Gran Telescopio Canarias

Context. Known Mars trojans could be primordial small bodies that have remained in their present-day orbits for the age of the Solar System. Their orbital distribution is strongly asymmetric; there are over a dozen objects at the L5 point and just one at L4, (121514) 1999 UJ7. Most L5 trojans appear to form a collision-induced asteroid cluster, known as the Eureka family. Asteroid 2023 FW14 was recently discovered and it has a robust orbit determination that may be consistent with a Mars trojan status. Aims. Our aim is determine the nature and dynamical properties of 2023 FW14. Methods. We carried out an observational study of 2023 FW14 to derive its spectral class using the OSIRIS camera spectrograph at the 10.4 m Gran Telescopio Canarias. We investigated its possible trojan resonance with Mars using direct N-body simulations. Results. The reflectance spectrum of 2023 FW14 is not compatible with the olivine-rich composition of the Eureka family; it also does not resemble the composition of the Moon, although (101429) 1998 VF31 does. The Eureka family and 101429 are at the L5 point. The spectrum of 2023 FW14 is also different from two out of the three spectra in the literature of the other known L4 trojan, 121514, which are of C-type. The visible spectrum of 2023 FW14 is consistent with that of an X-type asteroid, as is the third spectrum of 121514. Our calculations confirm that 2023 FW14 is the second known L4 Mars trojan although it is unlikely to be primordial; it may remain in its present-day “tadpole” path for several million years before transferring to a Mars-crossing orbit. It might be a fragment of 121514, but a capture scenario seems more likely. Conclusions. The discovery of 2023 FW14 suggests that regular Mars-crossing asteroids can be captured as temporary Mars trojans.

Direct N-body Simulations of Satellite Formation around Small Asteroids: Insights from DART’s Encounter with the Didymos System

We explore binary asteroid formation by spin-up and rotational disruption considering the NASA DART mission's encounter with the Didymos–Dimorphos binary, which was the first small binary visited by a spacecraft. Using a suite of N-body simulations, we follow the gravitational accumulation of a satellite from meter-sized particles following a mass-shedding event from a rapidly rotating primary. The satellite’s formation is chaotic, as it undergoes a series of collisions, mergers, and close gravitational encounters with other moonlets, leading to a wide range of outcomes in terms of the satellite's mass, shape, orbit, and rotation state. We find that a Dimorphos-like satellite can form rapidly, in a matter of days, following a realistic mass-shedding event in which only ∼2%–3% of the primary's mass is shed. Satellites can form in synchronous rotation due to their formation near the Roche limit. There is a strong preference for forming prolate (elongated) satellites, although some simulations result in oblate spheroids like Dimorphos. The distribution of simulated secondary shapes is broadly consistent with other binary systems measured through radar or lightcurves. Unless Dimorphos's shape is an outlier, and considering the observational bias against lightcurve-based determination of secondary elongations for oblate bodies, we suggest there could be a significant population of oblate secondaries. If these satellites initially form with elongated shapes, a yet-unidentified pathway is needed to explain how they become oblate. Finally, we show that this chaotic formation pathway occasionally forms asteroid pairs and stable triples, including coorbital satellites and satellites in mean-motion resonances.

Hot Jupiter formation in dense star clusters

ABSTRACT Hot Jupiters (HJ) are defined as Jupiter-mass exoplanets orbiting around their host star with an orbital period < 10 d. It is assumed that HJ do not form in-situ but ex-situ. Recent discoveries show that star clusters contribute to the formation of HJ. We present direct N-body simulations of planetary systems in star clusters and analyse the formation of HJ in them. We combine two direct N-body codes: nbody6++gpu for the dynamics of dense star clusters with 32 000 and 64 000 stellar members and lonelyplanets used to follow 200 identical planetary systems around solar mass stars in those star clusters. We use different sets with three, four, or five planets and with the innermost planet at a semimajor axis of 5 or 1 au and follow them for 100 Myr in our simulations. The results indicate that HJs are generated with high efficiency in dense star clusters if the innermost planet is already close to the host star at a semimajor axis of 1 au. If the innermost planet is initially beyond a semimajor axis of 5 au, the probability of a potential HJ ranges between 1.5 and 4.5 per cent. Very dense stellar neighbourhoods tend to eject planets rather than forming HJs. A correlation between HJ formation and angular momentum deficit is not witnessed. Young HJs (tage < 100 Myr) have only been found, in our simulations, in planetary systems with the innermost planet at a semimajor axis of 1 au.

ESTIMATING GAMMA-RAY FLUX FROM MILLISECOND PULSARS ORIGINATING IN GLOBULAR CLUSTERS NEAR THE GALACTIC CENTER

In this study, we investigate the contribution of millisecond pulsars (MSPs) to the gamma-ray excess observed in the Galactic Center by analyzing data from high-resolution direct N-body simulations of six globular clusters (GCs) that experience close encounters with the nuclear star cluster. Using the φ-GPU code, we tracked the orbits of individual neutron stars (NSs) formed during the simulations, assuming a fraction of these NSs evolve into MSPs. Our model includes state-of-the-art single stellar evolution code including prescription for neutron star formation. We estimated the gamma-ray flux from these MSPs, considering known values for their gamma-ray emission. Our results show that MSPs originating from the six modeled GCs contribute a small but non-negligible fraction of the observed gamma-ray flux. This finding suggests that the actual gamma-ray flux from MSPs could be much higher when considering the entire population of GCs, potentially significantly contributing to the gamma-ray excess. This study highlights the importance of considering MSPs in the Galactic Center, originating from nearby globular clusters, as a potential source of the observed gamma-ray excess. Future work will involve more sophisticated simulations incorporating binary stellar evolution and comparing the fraction of MSPs in observed GCs to refine our models and improve the accuracy of our estimates.

Enabling attractive-repulsive potentials in binary-collision-approximation monte-carlo codes for ion-surface interactions

Binary Collision Approximation (BCA) codes for ion-material interactions, such as SRIM, Tridyn, F-TRIDYN, and SDtrimSP, have historically been limited to screened Coulomb potentials even at low energies due to the difficulty in numerically solving the Distance of Closest Approach (DOCA) problem for attractive-repulsive potentials. Techniques such as direct n-body simulation or modifications to Newton’s method are either prohibitively costly or not guaranteed to work for all potentials. Advanced rootfinding techniques, such as companion matrix solvers, offer a solution. For many attractive-repulsive potentials, however, a companion matrix cannot be used directly, because there is no way to put the associated functions into a monomial basis form. A complementary technique is proxy rootfinding—by finding the best-fit polynomial approximant of a function, the zeros of the approximant can be guaranteed to be close to the zeros of the function. Using the Chebyshev basis and grid offers additional guarantees with regards to the quality of the approximation, the speed of convergence, and the avoidance of Runge’s phenomenon. By finding Chebyshev interpolants and using the Chebyshev-Frobenius companion matrix, the zeros of any real function on a bounded domain can be found. Here we show that using an Adaptive Chebyshev Proxy Rootfinder with Automatic Subdivision (ACPRAS) with appropriate scaling functions, numerical issues presented by attractive-repulsive potentials, including those of scale, can be handled. Using these techniques, we show that it is possible to include any physically reasonable interatomic potential in a BCA code, and to guarantee correctness of the resulting scattering angle calculations.

Exploring the Origin of Stars on Bound and Unbound Orbits Causing Tidal Disruption Events

Tidal disruption events (TDEs) provide a clue to the properties of a central supermassive black hole (SMBH) and an accretion disk around it, and to the stellar density and velocity distributions in the nuclear star cluster surrounding the SMBH. Deviations of TDE light curves from the standard occurring at a parabolic encounter with the SMBH depend on whether the stellar orbit is hyperbolic or eccentric and the penetration factor (β, the tidal disruption radius to the orbital pericenter ratio). We study the orbital parameters of bound and unbound stars being tidally disrupted by comparison of direct N-body simulation data with an analytical model. Starting from the classical steady-state Fokker–Planck model of Cohn & Kulsrud, we develop an analytical model of the number density distribution of those stars as a function of orbital eccentricity (e) and β. To do so, fittings of the density and velocity distribution of the nuclear star cluster and of the energy distribution of tidally disrupted stars are required and obtained from N-body data. We confirm that most of the stars causing TDEs in a spherical nuclear star cluster originate from the full loss-cone region of phase space, derive analytical boundaries in eccentricity-β space, and find them confirmed by N-body data. Since our limiting eccentricities are much smaller than critical eccentricities for full accretion or the full escape of stellar debris, we conclude that those stars are only very marginally eccentric or hyperbolic, close to parabolic.

Making hot Jupiters in stellar clusters – II. Efficient formation in binary systems

ABSTRACT Observations suggested that the occurrence rate of hot Jupiters (HJs) in open clusters is largely consistent with the field ($\sim 1{{\ \rm per\ cent}}$) but in the binary-rich cluster M67, the rate is $\sim 5{{\ \rm per\ cent}}$. How does the cluster environment boost HJ formation via the high-eccentricity tidal migration initiated by the extreme-amplitude von Zeipel–Lidov–Kozai (XZKL) mechanism forced by a companion star? Our analytical treatment shows that the cluster’s collective gravitational potential alters the companion’s orbit slowly, which may render the star–planet–companion configuration XZKL-favourable. We have also performed direct Gyr N-body simulations of the star cluster evolution and XZKL of planets’ orbit around member stars. We find that an initially single star may acquire a companion star via stellar scattering and the companion may enable XZKL in the planets’ orbit. Planets around an initially binary star may also be XZKL-activated by the companion. In both scenarios, the companion’s orbit has likely been significantly changed by stellar scattering and the cluster potential before XZKL occurs. Across different cluster models, 0.8–3 per cent of the planets orbiting initially single stars have experienced XZKL while the fraction is 2–26 per cent for initially binary stars. Around a star that is binary at 1 Gyr, 13–32 per cent of its planets have undergone XZKL, and combined with single stars, the overall XZKL fraction is 3–21 per cent, most affected by the cluster binarity. If 10 per cent of the stars in M67 host a giant planet, our model predicts an HJ occurrence rate of $\sim 1{{\ \rm per\ cent}}$. We suggest that HJ surveys target old, high-binarity, not-too-dense open clusters and prioritize wide binaries to maximize HJ yield.

Dynamical formation of Gaia BH1 in a young star cluster

ABSTRACT Gaia BH1, the first quiescent black hole (BH) detected from Gaia data, poses a challenge to most binary evolution models: its current mass ratio is ≈0.1, and its orbital period seems to be too long for a post-common envelope system and too short for a non-interacting binary system. Here, we explore the hypothesis that Gaia BH1 formed through dynamical interactions in a young star cluster (YSC). We study the properties of BH-main sequence (MS) binaries formed in YSCs with initial mass 3 × 102–3 × 104 M⊙ at solar metallicity, by means of 3.5 × 104 direct N-body simulations coupled with binary population synthesis. For comparison, we also run a sample of isolated binary stars with the same binary population synthesis code and initial conditions used in the dynamical models. We find that BH-MS systems that form via dynamical exchanges populate the region corresponding to the main orbital properties of Gaia BH1 (period, eccentricity, and masses). In contrast, none of our isolated binary systems match the orbital period and MS mass of Gaia BH1. Our best-matching Gaia BH1-like system forms via repeated dynamical exchanges and collisions involving the BH progenitor star, before it undergoes core collapse. YSCs are at least two orders of magnitude more efficient in forming Gaia BH1-like systems than isolated binary evolution.

Planetesimal Accretion at Short Orbital Periods

Abstract Formation models in which terrestrial bodies grow via the pairwise accretion of planetesimals have been reasonably successful at reproducing the general properties of the Solar System, including small-body populations. However, planetesimal accretion has not yet been fully explored in the context of the wide variety of recently discovered extrasolar planetary systems, in particular those that host short-period terrestrial planets. In this work, we use direct N-body simulations to explore and understand the growth of planetary embryos from planetesimals in disks extending down to ≃1 day orbital periods. We show that planetesimal accretion becomes nearly 100% efficient at short orbital periods, leading to embryo masses that are much larger than the classical isolation mass. For rocky bodies, the physical size of the object begins to occupy a significant fraction of its Hill sphere toward the inner edge of the disk. In this regime, most close encounters result in collisions, rather than scattering, and the system does not develop a bimodal population of dynamically hot planetesimals and dynamically cold oligarchs, as is seen in previous studies. The highly efficient accretion seen at short orbital periods implies that systems of tightly packed inner planets should be almost completely devoid of any residual small bodies. We demonstrate the robustness of our results to assumptions about the initial disk model, and we also investigate the effects that our simplified collision model has on the emergence of this non-oligarchic growth mode in a planet-forming disk.

Improving initialization and evolution accuracy of cosmological neutrino simulations

Neutrino mass constraints are a primary focus of current and future large-scale structure (LSS) surveys. Non-linear LSS models rely heavily on cosmological simulations — the impact of massive neutrinos should therefore be included in these simulations in a realistic, computationally tractable, and controlled manner. A recent proposal to reduce the related computational cost employs a symmetric neutrino momentum sampling strategy in the initial conditions. We implement a modified version of this strategy into the Hardware/Hybrid Accelerated Cosmology Code (HACC) and perform convergence tests on its internal parameters. We illustrate that this method can impart \U0001d4aa(1%) numerical artifacts on the total matter field on small scales, similar to previous findings, and present a method to remove these artifacts using Fourier-space filtering of the neutrino density field. Moreover, we show that the converged neutrino power spectrum does not follow linear theory predictions on relatively large scales at early times at the 15% level, prompting a more careful study of systematics in particle-based neutrino simulations. We also present an improved method for backscaling linear transfer functions for initial conditions in massive neutrino cosmologies that is based on achieving the same relative neutrino growth as computed with Boltzmann solvers. Our self-consistent backscaling method yields sub-percent accuracy in the total matter growth function. Comparisons for the non-linear power spectrum with the Mira-Titan emulator at a neutrino mass of mν = 0.15 eV are in very good agreement with the expected level of errors in the emulator and in the direct N-body simulation.

Dynamics of intermediate mass black holes in globular clusters

Context. We recently introduced a new method for simulating collisional gravitational N-body systems with approximately linear time scaling with N. Our method is based on the multi-particle collision (MPC) scheme, previously applied in fluid dynamics and plasma physics. We were able to simulate globular clusters with a realistic number of stellar particles (at least up to several times 106) on a standard workstation. Aims. We simulated clusters hosting an intermediate mass black hole (IMBH), probing a broad range of BH-cluster and BH–average-star mass ratios, unrestricted by the computational constraints that affect direct N-body codes. Methods. We set up a grid of hybrid particle-in-cell-MPC N-body simulations using our implementation of the MPC method, MPCDSS. We used either single mass models or models with a Salpeter mass function (a single power law with an exponent of −2.35), with the IMBH initially sitting at the centre. The force exerted by and on the IMBH was evaluated with a direct sum scheme with or without softening. For all simulations we measured the evolution of the Lagrangian radii and core density and velocity dispersion over time. In addition, we also measured the evolution of the velocity anisotropy profiles. Results. We find that models with an IMBH undergo core collapse at earlier times, the larger the IMBH mass the shallower they are, with an approximately constant central density at core collapse. The presence of an IMBH tends to lower the central velocity dispersion. These results hold independently of the mass function of the model. For the models with Salpeter MF, we observed that equipartition of kinetic energies is never achieved, even long after core collapse. Orbital anisotropy at large radii appears to be driven by energetic escapers on radial orbits, triggered by strong collisions with the IMBH in the core. We measured the wander radius, that is the distance of the IMBH from the centre of mass of the parent system over time, finding that its distribution has positive kurtosis. Conclusions. Among the results we obtained, which mostly confirm or extend previously known trends that had been established over the range of parameters accessible to direct N-body simulations, we underline that the leptokurtic nature of the IMBH wander radius distribution might lead to IMBHs presenting as off-centre more frequently than expected, with implications on observational IMBH detection.