EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 693 (January 2025) A&A, 693 (2025) A166 Full HTML
Open Access
Issue
A&A
Volume 693, January 2025
Article Number A166
Number of page(s) 10
Section Numerical methods and codes
DOI https://doi.org/10.1051/0004-6361/202452441
Published online 14 January 2025

© The Authors 2025

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

1 Introduction

The vast majority of young stars are found in or near star-forming regions that contain up to hundreds of thousands of stars (e.g., Lada & Lada 2003). Most of these star-forming regions disperse within several tens of millions of years, while others remain bound and become the open clusters and globular clusters that are seen in the Galaxy today (e.g., de Grijs & Parmentier 2007). Isotope studies of meteorites suggest that even our own Solar System was once part of a star cluster that has long dispersed (Looney et al. 2006; Portegies Zwart 2009). Given that most stars are formed in crowded stellar environments (e.g., Parker 2020), investigating the evolution of planetary systems and their constituents in such environments is worthwhile. Stellar encounters with planetary systems, as well as gravitational interactions between the constituents of a planetary system, can result in the ejection of planets, moons, asteroids, and comets, resulting in free-floating planetary debris in a star cluster. Hereafter, we collectively refer to bodies with masses far below those of stars as massless particles (MLPs).

Upon ejection from their planetary system, MLPs typically obtain a velocity at infinity in the range 0.1–10 km s−1. The velocity at infinity is higher for prompt ejections (as a direct result of star–planet encounters) than for delayed ejection (as a result of a planet–planet scattering; see Malmberg et al. 2011). When an MLP is ejected with a speed above the star cluster’s local escape velocity, it will escape from the star cluster. When it is ejected from its host planetary system with a speed below the star cluster’s local escape velocity, it remains part of the star cluster and becomes an MLP (e.g., Wu et al. 2023, 2024). We show in our study that the trajectories of such MLPs generally follow the gravitational potential of the star cluster. In specific cases, a free-floating MLP can also be captured by another star (e.g., Kouwenhoven et al. 2010; Malmberg et al. 2011; Moeckel & Clarke 2011; Parker & Quanz 2012; Perets & Kouwenhoven 2012; Portegies Zwart et al. 2021). Such captured bodies may be common (e.g., Siraj & Loeb 2020) and can potentially be identified through their orbital parameter distributions (e.g., Siraj & Loeb 2019).

Directly observing MLPs is extremely challenging, but extrapolations and surveys have suggested that MLPs are abundant in the Galactic disc. Initial estimation ranged from two Jupiter-mass MLPs for every main-sequence star (Sumi et al. 2011) to 105 MLPs with masses between 10−8 M⊙ and 10−2 M⊙ per star (Strigari et al. 2012). Most observed MLPs are detected with microlensing (Mao & Paczynski 1991; Gould & Loeb 1992; Abe et al. 2004; Beaulieu et al. 2006; Gaudi 2012). However, only a few MLPs have been unambiguously detected (e.g. CFBDSIR 21490403 and PSO J318-22; see Delorme et al. 2012; Liu et al. 2013, respectively). Most of the other MLPs have uncertain mass estimates, which prevents a confirmation. As of today, 26 observed candidates have been detected. The comprehensive analysis by Mróz et al. (2017) finds that Sumi et al. (2011) greatly overestimates the abundance of free-floating Jupiter-mass planets in the Galactic field and points out that the findings of Sumi et al. (2011) are inconsistent with the current planet formation theories (Veras & Raymond 2012; Ma et al. 2016). The former findings are also inconsistent with surveys of young clusters (Peña Ramírez et al. 2016; Ramirez & Kaltenegger 2017; Scholz et al. 2012; Mužić et al. 2015). Using a much larger sample of microlensing events, Mróz et al. (2017) revealed a significantly lower abundance of roughly one free-floating Jupiter-mass planet or wide-orbit Jupiter-mass planet for every four mainsequence stars in the Galactic field. More recently, a survey on stellar associations by Miret-Roig (2023) using the technique of microlensing found thousands of free-floating planet candidates.

In this work, we present the code NBODY6++GPU-MASSLESS, which can efficiently model the dynamical evolution of large populations of MLPs in star clusters. In this first application of the code, we consider the dynamical evolution of a population of free-floating MLPs in star clusters. We also investigate the kinematics of MLPs and their relation to the mass segregation process, changing the number density of stars in our star cluster models. Our findings are applicable to any population of bodies that can be approximated as massless, including comets, asteroids, planetesimals, and free-floating planets.

This paper is organised as follows. In Sect. 2 we describe the code NBODY6++GPU-MASSLESS, our numerical method, and initial conditions. We discuss the numerical performance of our code, in comparison with NBODY6++GPU in Sect. 3. We then discuss the dynamical evolution of the star clusters and the freefloating MLP populations in Sect. 4. Finally, we summarise our results and discuss the implications of our findings in Sect. 5.

2 NBODY6++GPU-MASSLESS

2.1 Theoretical background

In addition to containing to stars, star clusters also hold a vast number of smaller bodies, such as planets, moons, asteroids, comets, and planetesimals. These particles are so abundant that it is impractical to directly include them as particles in N-body simulations. Moreover, since these particles have a negligible mass, it is possible to neglect their gravitational force on the other bodies in the system. MLPs experience the gravitational forces of the massive bodies in a system, but they do not exert gravitational forces on any of the other bodies. NBODY6++GPU-MASSLESS makes use of this property of MLPs, which allows for a significant increase in the speed of simulations.

The wall-clock time of an N-body simulation is dominated by the regular and irregular force calculations (Ahmad & Cohen 1973). If we consider a system of N = Ns + Np particles, where Ns is the number of massive particles, such as stars, and Np is the number of MLPs. The number of regular force calculations for each time step is then NF1∝N(N−1)/2≈N2/2,$\[N_{F 1} \propto N(N-1) / 2 \approx N^{2} / 2,\]$(1)

when N ≫ 1. With the special treatment of MLPs, the number of irregular force calculations for each regular time step is NF2∝Ns(Ns−1)/2+NpNs≈(Ns/2+Np)Ns,$\[N_{F 2} \propto N_{s}\left(N_{s}-1\right) / 2+N_{p} N_{s} \approx\left(N_{s} / 2+N_{p}\right) N_{s},\]$(2)

when Ns, Np ≫ 1.

This approach can lead to a significant speed-up when Np > Ns. When Np = αNs, then NF1/NF2 ≈ (1 + α)2/(1 + 2 α). When α = 0, then NF1/NF2 = 0; when α = 1, then NF1/NF2 ≈ 1.3; when α = 102, then NF1/NF2 ≈ 51; and when α = 104, then NF1/NF2 ≈ 5000. By distinguishing between massive particles and MLPs in the N-body code, it is possible to efficiently integrate star clusters with a large number of small bodies, provided that the assumption of negligible mass is not violated.

2.2 Implementation in NBODY6++GPU-MASSLESS

We present NBODY6++GPU-MASSLESS, a modified version of NBODY6++GPU (see Aarseth 1999, 2003; Spurzem 1999; Wang et al. 2015b, 2016; Kamlah et al. 2022; Spurzem & Kamlah 2023). In this updated version, we have implemented a fast and accurate treatment of large numbers of MLPs. Each particle in the code is assigned a flag that indicates whether it is treated as a massive particle or as an MLP.

When calculating the gravitational interaction between particles in NBODY6++GPU-MASSLESS, the gravitational forces exerted by MLPs are ignored. MLPs are therefore never included in the force loops of any other particle. MLPs are excluded from the regular and irregular force calculations for the stars. MLPs are also excluded from being the primary component of a Kustaanheimo-Stiefel (KS) regularisation (Kustaanheimo & Stiefel 1965) event, although they can participate as a companion. For details about the algorithms mentioned above, we refer to Aarseth (2003).

For each star in a selected system, we identify its nearest neighbour using the Ahmad & Cohen (1973) scheme. The candidates for KS regularisation have to satisfy two additional conditions in NBODY6++GPU and NBODY6++GPU-MASSLESS: (i) R · V < 0.1 GMR$\[\sqrt{G M R}\]$ and (ii) |apert|R2/GM < 0.25. Here, R is the scalar distance, R is the relative position, and V is the relative velocity. The term apert is the vectorial differential acceleration exerted by other perturbing particles, M is the sum of the masses of the two candidates, and G is the gravitational constant. For further details, we refer to Spurzem & Kamlah (2023).

From the masses, positions, and velocities of the particles, we subsequently derive the orbital parameters of each binary system. All receding stars and MLPs located beyond twice the star cluster’s tidal radius are treated as escapers (see, e.g., Aarseth 2003).

Stellar evolution and binary evolution are implemented following the prescriptions of the Cambridge stellar evolution package and their improvements (Eggleton et al. 1989; Hurley et al. 2000, 2002, 2005; Belczynski et al. 2007; Kamlah et al. 2022). We adopt level C, as described in Table A1 of Kamlah et al. (2022). For relativistic kicks and mergers, we use the methodology described in Arca Sedda et al. (2023, 2024a,b).

3 Numerical performance

The wall-clock time consumption for a system consisting of N gravitationally interacting particles scales approximately as ΔTCPU≈k1N2+k2NN¯n+k3+k4N2+ΔTKS+ΔTcomm.$\[\Delta T_{\mathrm{CPU}} \approx k_{1} N^{2}+k_{2} N \bar{N}_{n}+k_{3}+k_{4} N^{2}+\Delta T_{\mathrm{KS}}+\Delta T_{\text {comm}}.\]$(3)

Here, k1 denotes the regular force computation coefficient, k2 denotes the irregular (neighbour) force coefficient, and N¯n$\[\bar{N}_{n}\]$ is the average neighbour number, which is typically 50–200 (Ahmad & Cohen 1973). The coefficient k3 represents the computational time that is (mostly) independent of N and N¯n$\[\bar{N}_{n}\]$ (for example, processes such as simulation parameter updates; see Huang et al. 2016, for details). The time required for diagnostics (e.g. the energy conservation checks) is governed by the coefficient k4. The term ΔTKS is the time required for integrating two-body KS regularisation pairs (e.g., Kustaanheimo & Stiefel 1965), few-body chains, and hierarchical systems. Finally, ΔTcomm represents the overhead term that accounts for parallel computation.

The contribution of the coefficients k3 and k4 to the wallclock time is negligible for our simulations. As we carry out energy checks every few million years, k4 ≪ k1. Moreover, ΔTKS is small (see also the KS implementation of MLPs above), and ΔTcomm is negligible, so for a star cluster with N = Ns + Np bodies, Eq. (3) thus reduces to ΔTCPU≈k1N2+k2NN¯n=k1(Ns2+2NsNp+Np2)+k2(Ns+Np)N¯n.$\[\Delta T_{\mathrm{CPU}} \approx k_{1} N^{2}+k_{2} N \bar{N}_{n}=k_{1}\left(N_{\mathrm{s}}^{2}+2 N_{\mathrm{s}} N_{\mathrm{p}}+N_{\mathrm{p}}^{2}\right)+k_{2}\left(N_{\mathrm{s}}+N_{\mathrm{p}}\right) \bar{N}_{n}.\]$(4)

As we disable the interaction between MLPs, the Np2$\[N_{p}^{2}\]$ and Ns2$\[N_{s}^{2}\]$ terms in the force calculations can be neglected. Also, the gravitational influence from the MLPs on the massive particles is disabled, so the term NsNp is reduced by half. With these modifications, the above expression then reduces to ΔTCPU≈(k1Ns+k2N¯n)Np+k1Ns2+k2NsN¯n.$\[\Delta T_{\mathrm{CPU}} \approx\left(k_{1} N_{\mathrm{s}}+k_{2} \bar{N}_{n}\right) N_{\mathrm{p}}+k_{1} N_{\mathrm{s}}^{2}+k_{2} N_{\mathrm{s}} \bar{N}_{n}.\]$(5)

Our test simulations were carried out with star clusters initialised with identical parameters except for the number of stars and MLPs. Each simulation was evolved for ten N-body units. We first modelled a star cluster without MLPs. We then added 10% of the total number of stars for each model until the number of MLPs was equal to the number of stars. We used one set of models with 64 000 stars and one set of models with 128 000 stars. We carried out our simulations on a workstation with an Intel CPU (i9–9960X CPU 3.00 GHz) and two NVIDIA GPU cards (GeForce RTX 2080 Ti). The wall-clock simulation time, ΔTCPU, for different values of Ns and Np is shown in Fig. 1 from staronly simulations up to Ns = Np. Figure 2 shows the distribution of the wall-clock time consumption over different components for the model with 64 000 stars and 64 000 MLPs and the model with 128 000 stars. The distribution on the pie chart is relative to the total time of the simulation and not to the seconds used in the simulation. The regular force calculations, irregular force calculations, and the adjustment procedures provide the largest contributions to the wall-clock time consumption.

The main difference in the computational time of the two simulations is related to the regular and irregular forces, as the other processes are also related to hardware properties and auxiliary softwares used on the machine. The regular forces take 13% less time in the star-MLP simulations. These forces are updated less frequently than the irregular forces, so the difference between the two simulations becomes more noticeable with a larger number of N-body temporal units. The main contribution resides in the irregular forces, where the star and MLP simulation takes only ≈30% of the time compared to the star-only simulation. The entire simulation, which lasts for ten N-body units, has a net gain in total computational time of ≈600 s (see Fig. 1). Our comparison between NBODY6++GPU and NBODY6++GPU-MASSLESS indicates that the codes provide the same performance (apart from statistical fluctuations) for star-only clusters, but the special treatment of MLPs in NBODY6++GPU-MASSLESS significantly speeds up the integration (depending on the number of MLPs relative to the number of stars). As planets, comets, and other debris can be safely treated as MLPs, the total number of such particles and their kinematic properties do not affect the dynamical evolution of the stellar population. Also, since MLPs do not interact with each other, all of our results can be easily scaled up to star clusters with a much larger number of MLPs.

Finally, in Fig. 3, we present a comparison between the standard version of NBODY6++GPU and the code presented in this paper, NBODY6++GPU-MASSLESS. We carried out simulations for model C12.8k (see Sect. 4.1) for 300 Myrs. The simulations have identical initial conditions, and the simulations have the same initial conditions for all particles, using a Mercury mass for all the MLPs. The only difference is that in the standard NBODY6++GPU version, the MLPs are included in the calculations of the regular and irregular forces. The difference in the Lagrangian radii is minimal, which justifies our treatment of MLPs.

thumbnail Fig. 1

Wall-clock time (TCPU) consumption for carrying out star cluster simulations with different numbers of stars (Ns) and different numbers of MLPs (Np). Each calculation was carried out for ten N-body units. Simulations were performed for Ns = 64 000 (blue curve) and Ns = 128 000 (orange curve) and with an initial MLP-to-star ratio (Np/Ns) ranging from zero to unity (blue dots).
thumbnail Fig. 2

Wall-clock calculation times for different processes in the simulation with Ns = 64 000 and Np = 64 000 (top pie chart) and Ns = 128 000 and Np = 0 (bottom pie chart) in Fig. 1. The chart presents the wall-clock time for time intervals of the total simulation time (ten code units). Different numerical processes are indicated by colour: regular force calculations (blue), irregular force calculations (orange), adjustment processes (green), and others (red). It is immediately noticeable that the irregular force calculations consistently take less time overall in the simulation with massless particles.
thumbnail Fig. 3

Comparison between the Lagrangian radii of the MLPs in the C12.8k model using the standard version of NBODY6++GPU (in orange) and using NBODY6++GPU-MASSLESS (in blue). The C12.8k model is also used for the main simulations (see next sections). The curves show the 0.1 (bottom), 1, 10, 50, and 90% (top) Lagrangian radii for both models. The difference between the two models is minimal and is a consequence of the stochastic nature of the clusters and numerical noise.

4 Dynamical evolution of MLP populations

4.1 Initial conditions

We numerically studied the evolution of free-floating MLPs in star clusters by modelling three sets of star cluster simulations, which we hereafter refer to as models C12.8k, C64.8k, and C128k. The choice of this particular set of models allowed us to determine how (i) the different stellar densities affect the dynamical evolution of MLPs and (ii) which dynamical processes affect the MLPs and to what degree they are affected.

The initial conditions for each of the three models are summarised in Table 1. These models represent star clusters containing Ns = 12 800, Ns = 64 800, and Ns = 128 000 stars, respectively. The stellar positions and velocities were drawn from a Plummer (1911) distribution in virial equilibrium. All models were initialised with virial radii of rvir = 1 pc, corresponding to an initial half-mass radii of rhm ≈ 0.77 pc. Stellar masses are drawn from the Kroupa (2001) initial mass function in the mass range 0.08–150 M⊙. The total initial masses of the three clusters are Mcl ≈ 7.45 × 103 M⊙, 3.7 × 104 M⊙, and 7.45 × 104 M⊙. The models do not include primordial binaries, and we do not include primordial mass segregation. The initial velocity distributions of the stellar and MLP populations are isotropic. We note that an initially anisotropic velocity distribution strongly affects mass segregation and relaxation in star clusters (e.g., Pavlík & Vesperini 2022; Tiongco et al. 2022; Livernois et al. 2022). The star clusters evolve in an external tidal field corresponding to that of a star cluster in a solar orbit in the Milky Way galaxy. The Milky Way was modelled as a point mass of MG = 9.56 × 1010 M⊙, with the star cluster in a circular orbit with radius RG = 8.5 kpc. The corresponding tidal radius (rt) of a cluster can then be estimated using rt ≈ (Mcl/3MG)1/3 RG.

In addition to the stars, each model containing N stars also contained N MLPs. The MLPs were assigned positions and velocities with distributions that are statistically identical to those of the stars. We note that in realistic star clusters, ejections of MLPs from their host planetary system can initially lead to a slightly super-virial population of MLPs in the star cluster. However, this population rapidly virialises after a brief epoch of rapid escape at early times; we refer the reader to Wang et al. (2015a) for a discussion on this issue. As the MLPs do not exert a gravitational force on the other bodies, they do not affect the evolution of the stellar population and the other MLPs (see Sect. 2.2).

We evolved all models for 300 Myr, which corresponds to roughly 5trh for model C128k. Model C64.8k is initially around five times denser than model C12.8k, and model C128k is initially about ten times denser than model C12.8k.

4.2 Timescales of star cluster evolution

The global evolution of a star cluster in virial equilibrium is governed by its initial crossing time, tcr; its initial half-mass relaxation time, trlx; and its mass segregation timescale, tms. The initial crossing times for the three models are 0.11, 0.05, and 0.03 Myr, respectively (see Table 1). The half-mass two-body relaxation time (Spitzer 1987) is defined as trlx=0.138Nsln⁡Λ(rhm3GMcl)1/2.$\[t_{\mathrm{rlx}}=\frac{0.138 N_{\mathrm{s}}}{\ln \Lambda}\left(\frac{r_{\mathrm{hm}}^{3}}{G M_{\mathrm{cl}}}\right)^{1 / 2}.\]$(6)

Here, rhm is the half-mass radius, ln Λ ≈ ln(0.11Ns) is the Coulomb logarithm (Giersz & Spurzem 1994), and Ns is the number of stars. The mass segregation timescale can be expressed as tms=m¯mmaxtrlx,$\[t_{\mathrm{ms}}=\frac{\bar{m}}{m_{\max }} t_{\mathrm{rlx}},\]$(7)

where m¯$\[\bar{m}\]$ is the average stellar mass and mmax ≈ 96 M⊙ is the mass of the most massive star. The mass segregation timescale is roughly tms ≈ 0.006trlx, resulting in 0.165 Myr, 0.31 Myr, and 0.40 Myr for model C12.8k, C64.8k, and C128k, respectively. We note that the properties of the MLP population do not impact tcr, trlx, and tms, as MLPs do not affect the dynamics of the stellar population.

Table 1

Initial conditions of the star cluster models.

4.3 Lagrangian radii and mass segregation

Lagrangian radii provide a powerful tool for analysing the global evolution of a star cluster, as well as for characterising mass segregation. The evolution of the Lagrangian radii for the three models is shown in Fig. 4 for Lagrangian radii ranging from the 0.1% mass shell to the 90% mass shell.

The Lagrangian radii for the stellar population were calculated using the total cluster mass at each time. The Lagrangian radii for the population of MLPs were calculated from the number of MLPs enclosed in the Lagrangian shell at that time (instead of the mass, since MLPs are massless).

During the core collapse phase of a star cluster, strong gravitational scattering events result in the ejection of stars to (and beyond) the outskirts of the star cluster. Simultaneously, stellar evolution reduces the total mass of the star cluster. As a result of stellar ejections and stellar mass loss, the gravitational potential of the star cluster gradually decreases, and MLPs slowly migrate towards the outskirts of the star cluster. This results in a divergence between the Lagrangian radii for the stars and the Lagrangian radii for the MLPs. As the expansion of the MLP population follows the gravitational potential of the star cluster, the Lagrangian radii for the MLPs are typically larger than the corresponding Lagrangian radii of the stellar population, except in the outskirts of the cluster, which is dominated by stars ejected from the core, primarily during the core collapse phase.

The evolution of the Lagrangian radii is similar in all models, but there are notable variations due to the differences in stellar number density of the star clusters. The inner regions display a significant change between the dynamical evolution of stars and MLPs. The innermost regions in Fig. 4 show a less pronounced core collapse for the Lagrangian radii of the MLPs when compared to the Lagrangian radii of the stars. The expansion of the population of MLPs in the inner regions is less pronounced when the stellar density is higher. In model C128k, for example, the 10% shell of the MLPs and stars is similar.

The 50% Lagrangian radii of the stars and MLPs evolve similarly but show a divergence after the core collapse. The stellar 50% Lagrangian radius grows faster than the corresponding MLP Lagrangian radius in model C12.8k.

The stellar Lagrangian 90% shell is particularly interesting. In dense models, after the initial core collapse, there is a larger ejection of mass in the cluster. The former leads to a significant expansion of the external regions due to the ejection of highenergy stars, particularly black holes and neutron stars, which is due to the stellar evolution in the first million years. The large expansion of the shell can be explained by the relatively mild collapse of the core in the inner regions. This causes high-energy stars to be ejected, leading to an earlier cooling of the core. The 90% Lagrangian shells of the MLPs follow a similar evolution in all models, regardless of the density of the star cluster. In particular, the 50% Lagrangian shell of the stars in model C128k becomes similar to the 90% Lagrangian shell of the MLPs, indicating that MLPs are more likely to remain bounded to the star cluster. The spike in the 90% shell of the MLPs Lagrangian radius in the middle panel of Fig. 4 is caused by an MLP that is ejected from the star cluster at the next time step.

More massive inner regions are more likely to bound the MLPs. Thus, they migrate much more slowly to the outer regions than the stars. The behaviour of the 90% shell is related to (i) the ejection of a large number of hard binaries, which make a large mass contribution to outer shells, and (ii) black holes and neutron stars, which also contribute to a less strong core collapse in the denser models.

In stellar systems containing two populations that have considerable differences in particle masses, the high-mass component tends to decouple from the low-mass component within a short period of time (see Khalisi et al. 2007). As the MLPs are massless, they do not participate in the energy equipartition process that is responsible for mass segregation of the stellar population. Free-floating MLPs thus have relatively unperturbed trajectories, and they generally follow the overall gravitational potential of the star cluster.

Close encounters between stars and MLPs do occur (see, e.g., Wang et al. 2015b). During a two-body scattering event between a star and an MLP, the motion of the star remains unperturbed. In the absence of other massive perturbers nearby, the MLP’s trajectory in the two-body encounter can be approximated with a hyperbolic orbit, during which the direction of its velocity vector changes, while the magnitude of the velocity vector remains constant (when evaluated at times sufficiently long before and after the encounter). The latter interaction therefore generally does not result in the ejection of an MLP, as its speed remains below the local escape velocity.

Close encounters involving three stars can result in the exchange of a substantial amount of kinetic energy and may lead to the ejection of one of the stars, usually the lowestmass star. Three-body encounters involving two stars and one MLP occur less frequently, as the trajectories are less affected by gravitational focusing (with respect to an encounter between three stars). When the three bodies approach each other at sufficiently close distances, this process contributes to the ejection of MLPs following gravitational slingshot events (see, e.g., Saslaw et al. 1974; Wang et al. 2015b). Three-body encounters involving one star and two MLPs can be treated as two separate two-body encounters, as the two MLPs do not interact with each other and therefore do not lead to ejection of any of the bodies from the system.

thumbnail Fig. 4

Evolution of the Lagrangian radii containing 0.1 (bottom), 1, 10, 50, and 90% (top) of the cluster mass for stars (in green) and MLPs (in red) for models C12.8k (top panel), C64.8k (middle panel), and C128k (bottom panel). The Lagrangian radii at each time are calculated using the total cluster mass of the star cluster at that time.
thumbnail Fig. 5

Evolution of the average stellar mass within the 0.1 (orange curve), 1 (red curve), 10 (green curve), 50 (blue curve), and 90% (black curve) Lagrangian shells for models C12.8k (top), C64.8k (middle), and C128k (bottom). The evolution in the first million years is dominated by stellar evolution and mass segregation, and it is steadier in all shells around 20 Myrs. The average mass of the inner Lagrangian shells changes more frequently due to the lower number of stars in the shells. This is more evident in the low density model.

4.4 Mass segregation

The presence of a mass spectrum in the stellar population results in rapid mass segregation (see Fig. 5). This process is a consequence of the tendency of the stellar population to achieve local energy equipartition, although the stellar population is unable to achieve energy equipartition through this process (see Parker et al. 2016). For details on the mass segregation process itself, we refer to the studies of Khalisi et al. (2007) and Allison et al. (2009).

Figure 5 shows the average stellar mass contained within several Lagrangian shells. The scatter in the curves representing the innermost Lagrangian radii is affected by small-number statistics. During the first million years, the tendency of stars to achieve local energy equipartition causes the most massive stars to sink to the centre, where they experience close encounters (e.g., Khalisi et al. 2007). This causes rapid variations in the average mass in the core. These strong fluctuations are gradually mediated by stellar evolution and expansion of the core. Throughout most of the simulations, the average mass in the core region is typically a factor of ten higher than that in the outskirts of the star cluster.

The average stellar mass within the 1% Lagrangian radius (blue curves in Fig. 5) increases substantially as the cluster evolves. The initial increase (t < 5 Myr) is due to mass segregation, and the substantial decrease at later times is caused by stellar evolution. Inspection of the output files of the simulations indicates that the central region is initially dominated by several massive stars. As most massive stars have evolved at times t ≳ 10 Myr, stellar evolution becomes less prominent, and the average mass in the core is dominated by stellar mass segregation. In the denser star cluster models, mass segregation occurs at later times, mass segregation is more prominent, and the outer regions abruptly stabilise at the time of the core collapse.

thumbnail Fig. 6

Distances and speeds of MLPs (left column) and stars (right column) at t = 0 Myr (top) and at t = 300 Myr (bottom) for model C128k. Colours indicate the number density of the data in the plot, as shown in the colourbar. For the methodology used, we refer to Flammini Dotti et al. (2023).

4.5 Escaping stars and MLPs

Figure 6 shows the radial distances and speeds of all stars and MLPs in the C128k model at times t = 0 Myr and t = 300 Myr. In the figure, stars are ejected from the core region of each cluster as a result of close encounters with other stars, and mass is lost due to stellar evolution. Most of the stars in the outskirts of the star cluster migrate to outer regions as a consequence of scattering with other stars and to a lesser extent due to the reduced gravitational potential of the star cluster following stellar evolution and stellar escapers. Most of the MLPs that migrate to the outskirts do so because of the decreasing gravitational potential of the star cluster and tend to escape over time, as the tidal radius of the star cluster shrinks, while the star cluster experiences mass loss due to stellar evolution and stellar escapers.

Figure 7 shows the cumulative number of escaping stars and MLPs over time. The distribution of escape speeds over time is shown in Fig. 8, and the evolution of the total cluster mass is shown in Fig. 9.

The escape rate of stars is higher during the initial phases of evolution due to the initial core collapse, and it gradually decreases as the cluster expands. The effect is definitely more prominent and lasts longer in dense star clusters. The subsequent expansion of the core, combined with the total mass, results in a longer relaxation time, implying that the two-body relaxation in the core, and consequently the stellar ejection rate, is suppressed (see Kouwenhoven et al. 2014). Therefore, stars are ejected mostly due to the core collapse at the start of the simulation and due to subsequent mass segregation. During the first 50 Myr, the ejection of high-energy stars (Fig. 8) is more prominent in denser star clusters. Most are neutron stars and black holes, which also indicates the prominent role of stellar evolution in the evolution of the star cluster.

The average speed of the ejected stars is 43.91, 107.14, and 134.56 km s−1 in models C12.8k, C64.8k, and C128k, respectively. The average speed of the ejected MLPs is 2.83, 6.52, and 9.42 km s−1 for models C12.8k, C64.8k, and C128k, respectively. The majority of the stellar escape events occur during the initial phases of the star cluster evolution. During this period, there are a large number of high-velocity escapers, which significantly impacts the average escape speed. In denser models, this region becomes even more populated. If we exclude these black holes and neutron stars from the average stellar escape speed, we obtain average escape speeds of 3.86, 8.26, and 11.56 km s−1 for models C12.8k, C64.8k, and C128k, respectively. The escaping MLPs thus have similar but slightly lower average escape speeds compared to the stellar escapers.

We did not find escaping MLPs with speeds greater than 102 km s−1. The majority of stars with such a high velocity are either binary systems or compact objects and are therefore unlikely to be found for an MLP. The zero mass of an MLP reflects in the total energy of the particle, which can only be enhanced to these velocities by interactions with stars or, more effectively, compact objects. To summarise, stars acquire large velocities also as a consequence of mass segregation and stellar evolution. Planets, however, acquire large velocities mostly through interactions with stars. Figure 8 shows that the escape speed distribution is mostly similar for stars and MLPs, with a somewhat broader escape speed distribution for the MLPs. After the first 50 Myr, the number of ejected stars is roughly 37, 42, and 50% of the total number of ejected stars at the end of the simulation. Therefore, these stars have a larger statistical weight on the average velocity. The stars are therefore mainly ejected from the star cluster in the first 50 Myrs due to the initial core collapse and stellar evolution. After the first 50 Myrs, the ejection of stars from the star cluster is due to either evaporation or ejection. The MLPs, instead, are ejected from the star cluster by either evaporation or ejection, and they are not related to the mechanism of mass segregation. This is another important result, as it shows there is no mass segregation impact on planet-mass objects.

Finally, we show the evolution of the total mass of each cluster in Fig. 9. During the 300 Myr of evolution, cluster models C12.8k, C64.8k, and C128k lose 25.81, 13.16, and 12.25% of their initial mass, respectively. The majority of mass is lost before 50 Myr for all models. Thus, stellar evolution greatly impacts the mass loss in denser star cluster models, which strongly reduces the mass loss. The mass loss is roughly linear with evaporation from the star cluster and high-energy stars for the remainder of the simulations.

To summarise, the ejection of high-energy stars impacts the loss of mass, while this contribution becomes less important for the remainder of the simulation. Moreover, denser star cluster models retain more mass due to their deeper gravitational potential wells.

thumbnail Fig. 7

The cumulative distribution of ejection times for MLPs (blue) and stars (orange) for models C12.8k (top), C64.8k (centre), and C128k (bottom).
thumbnail Fig. 8

The ejection velocity distribution of MLPs (blue) and stars (red) for models C12.8k (top), C64.8k (centre), and C128k (bottom).
thumbnail Fig. 9

Evolution of the total mass of the star cluster for all models. The total mass of the cluster gradually decreases over time due to stellar evolution and due to escaping stars. We note that the total mass is independent of the number of MLPs.

5 Conclusions

We have presented the first work done with the newly developed code NBODY6++GPU-MASSLESS, which is based on the direct N-body code NBODY6++GPU. NBODY6++GPU-MASSLESS is optimised for integrating star clusters with a large number of massless particles (MLPs), such as planets, comets, asteroids, and planetesimals. Using theoretical arguments and numerical performance tests, we have demonstrated that NBODY6++GPU-MASSLESS outperforms NBODY6++GPU significantly when a large number of MLPs are present. Through simulations of different star cluster models with a large number of MLPs, we analysed the dynamical evolution of the stellar and MLP populations, with an emphasis on the evolution of the Lagrangian radii, mass segregation, and the properties of escaping particles. Our main results are summarised as follows:

  1. We have developed the code NBODY6++GPU-MASSLESS, a modified version of NBODY6++GPU that allows for fast evolution of star clusters that contain a large number of MLPs. This was achieved through excluding the gravitational forces of MLPs on any of the other particles. We show that even when the number of MLPs is comparable to the number of stars, the wall-clock time is only moderately longer compared to the simulation without MLPs.

  2. Massless particles do not participate in the mass segregation process, as they are unable to change energy with the stellar population since they are test particles. In terms of their dynamical evolution, our study finds that (i) although the spatial distributions of the stellar population and the MLP population diverge over time, this process should not be confused with the classical mass segregation process for the stellar population. Unlike the stellar population, the kinematic evolution of a population of MLPs in a star cluster is primarily determined by the overall changes in the gravitational potential of the cluster, which is reduced over time due to stellar ejections and stellar evolution. This occurs when modelling MLPs as massless and when assigning them a small mass. Thus, modelling MLPs as massless is an adequate approximation (see Fig. 3).

  3. The dynamical evolution of the stellar population in star clusters with a higher initial stellar density is more strongly dominated by the initial core collapse, stellar evolution, and mass segregation, especially during the first 50 Myr. The population of MLPs, on the other hand, evolves similarly in all star cluster models.

  4. The cumulative number of escaping MLPs over time is smoother compared to that of the stellar escapers. There are two main reasons for this behaviour: (i) stellar evolution, which is especially represented by the high-energy stars ejected during the first 50 Myr (mainly neutron stars), and (ii) mass segregation, which inevitably causes the ejection of stars into the outer regions. This also tells us that those particles are mainly ejected through evaporation or ejection, but these processes are not influenced by mass segregation and stellar evolution.

  5. The average escape speeds of stars and MLPs are different. This is caused mainly by (i) stellar evolution, where high-velocity stars (such as neutron stars) are ejected, and (ii) density dependence, as the escape velocity is higher for denser star clusters. If we do not consider these high-energy stars during the first 50 Myrs (i.e. we ignore the contribution of stellar evolution to the stellar escape rate), then the ejection velocities of stars and MLPs are comparable. This also reflects what occurs during the Lagrangian evolution of the particles at the beginning of the simulation.

Though the parameter space covered by the models in our study is limited, our work provides useful insights into the differences between the dynamical evolution of the stellar population and the MLP population in a star cluster. We have not included primordial (stellar) binary systems in our simulations, nor have we included any planetary companions among star cluster members. These two features will be added in the subsequent papers in this series.

Acknowledgements

We are grateful to Long Wang for his support during the initial phase of this project. F.F.D. and R.S. acknowledge support by the German Science Foundation (DFG), priority program SPP 1992 “Exploring the diversity of extrasolar planets” (project Sp 345/22-1 and central visitor program). M.B.N.K. was supported by the National Natural Science Foundation of China (grant 11573004). F.F.D. and M.B.N.K. were supported by the Research Development Fund (grant RDF-16–01–16) of XJTLU. F.F.D. acknowledges support from the XJTLU postgraduate research scholarship. This publication was made possible through the support of a grant from the John Templeton Foundation and National Astronomical Observatories of the Chinese Academy of Sciences. P.B. thanks the support from the special program of the Polish Academy of Sciences and the U.S. National Academy of Sciences under the Long-term program to support Ukrainian research teams grant No. PAN.BFB.S.BWZ.329.022.2023. P.B. acknowledges the support by the National Science Foundation of China under grant NSFC No. 12473017.

References

  1. Aarseth, S. J. 1999, PASP, 111, 1333 [NASA ADS] [CrossRef] [Google Scholar]
  2. Aarseth, S. J. 2003, Gravitational N-Body Simulations, 1st edn. (Cambridge: Cambridge University Press), 430 [Google Scholar]
  3. Abe, F., Bennett, D. P., Bond, I. A., et al. 2004, Science, 305, 1264 [NASA ADS] [CrossRef] [Google Scholar]
  4. Ahmad, A., & Cohen, L. 1973, J. Comput. Physics, 12, 389 [NASA ADS] [CrossRef] [Google Scholar]
  5. Allison, R. J., Goodwin, S. P., Parker, R. J., et al. 2009, ApJ, 700, L99 [Google Scholar]
  6. Arca Sedda, M., Kamlah, A. W. H., Spurzem, R., et al. 2023, MNRAS, 526, 429 [NASA ADS] [CrossRef] [Google Scholar]
  7. Arca Sedda, M., Kamlah, A. W. H., Spurzem, R., et al. 2024a, MNRAS, 528, 5119 [NASA ADS] [CrossRef] [Google Scholar]
  8. Arca Sedda, M., Kamlah, A. W. H., Spurzem, R., et al. 2024b, MNRAS, 528, 5140 [CrossRef] [Google Scholar]
  9. Beaulieu, J.-P., Bennett, D. P., Fouqué, P., et al. 2006, Nature, 439, 437 [Google Scholar]
  10. Belczynski, K., Taam, R. E., Kalogera, V., Rasio, F. A., & Bulik, T. 2007, ApJ, 662, 504 [Google Scholar]
  11. de Grijs, R., & Parmentier, G. 2007, Chinese J. Astron. Astrophys., 7, 155 [NASA ADS] [CrossRef] [Google Scholar]
  12. Delorme, P., Gagné, J., Malo, L., et al. 2012, A&A, 548, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  13. Eggleton, P. P., Fitchett, M. J., & Tout, C. A. 1989, ApJ, 347, 998 [NASA ADS] [CrossRef] [Google Scholar]
  14. Flammini Dotti, F., Capuzzo-Dolcetta, R., & Kouwenhoven, M. B. N. 2023, MNRAS, 526, 1987 [NASA ADS] [CrossRef] [Google Scholar]
  15. Gaudi, B. S. 2012, ARA&A, 50, 411 [Google Scholar]
  16. Giersz, M., & Spurzem, R. 1994, MNRAS, 269, 241 [Google Scholar]
  17. Gould, A., & Loeb, A. 1992, ApJ, 396, 104 [Google Scholar]
  18. Huang, S.-Y., Spurzem, R., & Berczik, P. 2016, Res. Astron. Astrophys., 16, 11 [Google Scholar]
  19. Hurley, J. R., Pols, O. R., & Tout, C. A. 2000, MNRAS, 315, 543 [Google Scholar]
  20. Hurley, J. R., Tout, C. A., & Pols, O. R. 2002, MNRAS, 329, 897 [Google Scholar]
  21. Hurley, J. R., Pols, O. R., Aarseth, S. J., & Tout, C. A. 2005, MNRAS, 363, 293 [NASA ADS] [CrossRef] [Google Scholar]
  22. Kamlah, A., Leveque, A., Spurzem, R., et al. 2022, MNRAS, 511, 4060 [NASA ADS] [CrossRef] [Google Scholar]
  23. Khalisi, E., Amaro-Seoane, P., & Spurzem, R. 2007, MNRAS, 374, 703 [NASA ADS] [CrossRef] [Google Scholar]
  24. Kouwenhoven, M. B. N., Goodwin, S. P., Parker, R. J., et al. 2010, MNRAS, 404, 1835 [NASA ADS] [Google Scholar]
  25. Kouwenhoven, M. B. N., Goodwin, S. P., de Grijs, R., Rose, M., & Kim, S. S. 2014, MNRAS, 445, 2256 [NASA ADS] [CrossRef] [Google Scholar]
  26. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  27. Kustaanheimo, P., & Stiefel, E. 1965, J. Reine Angew. Math., 218, 204 [Google Scholar]
  28. Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57 [Google Scholar]
  29. Liu, H.-G., Zhang, H., & Zhou, J.-L. 2013, ApJ, 772, 13 [NASA ADS] [CrossRef] [Google Scholar]
  30. Livernois, A. R., Vesperini, E., Varri, A. L., et al. 2022, MNRAS, 512, 2584 [NASA ADS] [CrossRef] [Google Scholar]
  31. Looney, L. W., Tobin, J. J., & Fields, B. D. 2006, ApJ, 652, 1755 [NASA ADS] [CrossRef] [Google Scholar]
  32. Ma, S., Mao, S., Ida, S., Zhu, W., & Lin, D. N. C. 2016, MNRAS, 461, L107 [Google Scholar]
  33. Malmberg, D., Davies, M. B., & Heggie, D. C. 2011, MNRAS, 411, 859 [Google Scholar]
  34. Mao, S., & Paczynski, B. 1991, ApJ, 374, L37 [NASA ADS] [CrossRef] [Google Scholar]
  35. Miret-Roig, N. 2023, Ap&SS, 368, 17 [NASA ADS] [CrossRef] [Google Scholar]
  36. Moeckel, N., & Clarke, C. J. 2011, MNRAS, 415, 1179 [NASA ADS] [CrossRef] [Google Scholar]
  37. Mróz, P., Udalski, A., Skowron, J., et al. 2017, Nature, 548, 183 [Google Scholar]
  38. Mužić, K., Scholz, A., Geers, V. C., & Jayawardhana, R. 2015, ApJ, 810, 159 [Google Scholar]
  39. Parker, R. J. 2020, R. Soc. Open Sci., 7, 201271 [NASA ADS] [CrossRef] [Google Scholar]
  40. Parker, R. J., & Quanz, S. P. 2012, MNRAS, 419, 2448 [Google Scholar]
  41. Parker, R. J., Goodwin, S. P., Wright, N. J., Meyer, M. R., & Quanz, S. P. 2016, MNRAS, 459, L119 [NASA ADS] [CrossRef] [Google Scholar]
  42. Pavlík, V., & Vesperini, E. 2022, MNRAS, 515, 1830 [CrossRef] [Google Scholar]
  43. Peña Ramírez, K., Béjar, V. J. S., & Zapatero Osorio, M. R. 2016, A&A, 586, A157 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Perets, H. B., & Kouwenhoven, M. B. N. 2012, ApJ, 750, 83 [NASA ADS] [CrossRef] [Google Scholar]
  45. Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
  46. Portegies Zwart, S. F. 2009, ApJ, 696, L13 [NASA ADS] [CrossRef] [Google Scholar]
  47. Portegies Zwart, S., Torres, S., Cai, M. X., & Brown, A. G. A. 2021, A&A, 652, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Ramirez, R. M., & Kaltenegger, L. 2017, ApJ, 837, L4 [Google Scholar]
  49. Saslaw, W. C., Valtonen, M. J., & Aarseth, S. J. 1974, ApJ, 190, 253 [CrossRef] [Google Scholar]
  50. Scholz, A., Jayawardhana, R., Muzic, K., et al. 2012, ApJ, 756, 24 [NASA ADS] [CrossRef] [Google Scholar]
  51. Siraj, A., & Loeb, A. 2019, ApJ, 872, L10 [NASA ADS] [CrossRef] [Google Scholar]
  52. Siraj, A., & Loeb, A. 2020, arXiv e-prints [arXiv:2011.14900] [Google Scholar]
  53. Spitzer, L. 1987, Dynamical Evolution of Globular Clusters (Princeton: Princeton University Press) [Google Scholar]
  54. Spurzem, R. 1999, J. Comput. Appl. Math., 109, 407 [CrossRef] [Google Scholar]
  55. Spurzem, R., & Kamlah, A. 2023, Liv. Rev. Comput. Astrophys., 9, 3 [CrossRef] [Google Scholar]
  56. Strigari, L. E., Barnabè, M., Marshall, P. J., et al. 2012, MNRAS, 423, 1856 [NASA ADS] [CrossRef] [Google Scholar]
  57. Sumi, T., Kamiya, K., Bennett, D. P., et al. 2011, Nature, 473, 349 [Google Scholar]
  58. Tiongco, M. A., Vesperini, E., Varri, A. L., et al. 2022, MNRAS, 512, 1584 [NASA ADS] [CrossRef] [Google Scholar]
  59. Veras, D., & Raymond, S. N. 2012, MNRAS, 421, L117 [NASA ADS] [CrossRef] [Google Scholar]
  60. Wang, L., Kouwenhoven, M. B. N., Zheng, X., Church, R. P., & Davies, M. B. 2015a, MNRAS, 449, 3543 [NASA ADS] [CrossRef] [Google Scholar]
  61. Wang, L., Spurzem, R., Aarseth, S., et al. 2015b, MNRAS, 450, 4070 [NASA ADS] [CrossRef] [Google Scholar]
  62. Wang, L., Spurzem, R., Aarseth, S., et al. 2016, MNRAS, 458, 1450 [Google Scholar]
  63. Wu, K., Kouwenhoven, M. B. N., Spurzem, R., & Pang, X. 2023, MNRAS, 523, 4801 [NASA ADS] [CrossRef] [Google Scholar]
  64. Wu, K., Kouwenhoven, M. B. N., Dotti, F. F., & Spurzem, R. 2024, MNRAS, 533, 4485 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Initial conditions of the star cluster models.

In the text

All Figures

thumbnail Fig. 1

Wall-clock time (TCPU) consumption for carrying out star cluster simulations with different numbers of stars (Ns) and different numbers of MLPs (Np). Each calculation was carried out for ten N-body units. Simulations were performed for Ns = 64 000 (blue curve) and Ns = 128 000 (orange curve) and with an initial MLP-to-star ratio (Np/Ns) ranging from zero to unity (blue dots).
In the text
thumbnail Fig. 2

Wall-clock calculation times for different processes in the simulation with Ns = 64 000 and Np = 64 000 (top pie chart) and Ns = 128 000 and Np = 0 (bottom pie chart) in Fig. 1. The chart presents the wall-clock time for time intervals of the total simulation time (ten code units). Different numerical processes are indicated by colour: regular force calculations (blue), irregular force calculations (orange), adjustment processes (green), and others (red). It is immediately noticeable that the irregular force calculations consistently take less time overall in the simulation with massless particles.
In the text
thumbnail Fig. 3

Comparison between the Lagrangian radii of the MLPs in the C12.8k model using the standard version of NBODY6++GPU (in orange) and using NBODY6++GPU-MASSLESS (in blue). The C12.8k model is also used for the main simulations (see next sections). The curves show the 0.1 (bottom), 1, 10, 50, and 90% (top) Lagrangian radii for both models. The difference between the two models is minimal and is a consequence of the stochastic nature of the clusters and numerical noise.
In the text
thumbnail Fig. 4

Evolution of the Lagrangian radii containing 0.1 (bottom), 1, 10, 50, and 90% (top) of the cluster mass for stars (in green) and MLPs (in red) for models C12.8k (top panel), C64.8k (middle panel), and C128k (bottom panel). The Lagrangian radii at each time are calculated using the total cluster mass of the star cluster at that time.
In the text
thumbnail Fig. 5

Evolution of the average stellar mass within the 0.1 (orange curve), 1 (red curve), 10 (green curve), 50 (blue curve), and 90% (black curve) Lagrangian shells for models C12.8k (top), C64.8k (middle), and C128k (bottom). The evolution in the first million years is dominated by stellar evolution and mass segregation, and it is steadier in all shells around 20 Myrs. The average mass of the inner Lagrangian shells changes more frequently due to the lower number of stars in the shells. This is more evident in the low density model.
In the text
thumbnail Fig. 6

Distances and speeds of MLPs (left column) and stars (right column) at t = 0 Myr (top) and at t = 300 Myr (bottom) for model C128k. Colours indicate the number density of the data in the plot, as shown in the colourbar. For the methodology used, we refer to Flammini Dotti et al. (2023).
In the text
thumbnail Fig. 7

The cumulative distribution of ejection times for MLPs (blue) and stars (orange) for models C12.8k (top), C64.8k (centre), and C128k (bottom).
In the text
thumbnail Fig. 8

The ejection velocity distribution of MLPs (blue) and stars (red) for models C12.8k (top), C64.8k (centre), and C128k (bottom).
In the text
thumbnail Fig. 9

Evolution of the total mass of the star cluster for all models. The total mass of the cluster gradually decreases over time due to stellar evolution and due to escaping stars. We note that the total mass is independent of the number of MLPs.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.