EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 667 (November 2022) A&A, 667 (2022) A69 Full HTML
Open Access
Issue
A&A
Volume 667, November 2022
Article Number A69
Number of page(s) 9
Section Galactic structure, stellar clusters and populations
DOI https://doi.org/10.1051/0004-6361/202243579
Published online 07 November 2022

© P. Suin et al. 2022

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 and physical scenario

Understanding the evolution of star cluster has been the subject of many efforts over the past decades. Much work has been done to uncover how the internal effects of dynamical interactions (e.g., Spitzer & Shapley 1940; Chandrasekhar 1942; Lynden-Bell & Eggleton 1980; Spitzer 1987) and stellar evolution (e.g., Chernoff & Weinberg 1990; Portegies Zwart et al. 2007) impact its properties. Yet no cluster is truly isolated, and modelling the interaction with the environment, such as the large scale Galactic potential, spiral arms and high-density structure in the interstellar medium (ISM; e.g., Spitzer 1987; Terlevich 1987; Gnedin & Ostriker 1997; Gieles et al. 2006) also received great attention from the scientific community. Indeed, star clusters and associations originate in molecular clouds that are restructured through feedback from the stars, such as winds and explosive mass ejections (Lada & Lada 2003; Colín et al. 2013; Kruijssen et al. 2019; Adamo et al. 2020; Kim et al. 2021). Whatever the history of this interaction is, once free from the immediate natal gas, a young cluster still finds itself embedded in the larger scale, structurally complex environment from which the parent cloud emerged. This environment is dynamically chaotic since turbulence produces density and velocity fluctuations on a wide range of length scales (e.g., Vázquez-Semadeni 1994; Federrath & Klessen 2013; Chevance et al. 2020). These dense, evanescent turbulent structures generate a fluctuating potential that will harass the cluster and alter its evolution from that expected if it were isolated (Kruijssen et al. 2012). Before exiting from the formation region, the cluster has to pass through this denser environment. Here it will experience a variable potential that can distort its spatial distribution and increase its global velocity dispersion (see, e.g., Gieles et al. 2006).

In this work, we study how this gaseous environment affects a cluster’s dynamical evolution through numerical models that treat the hydrodynamics of the gas and the dynamics of the stars simultaneously. We focused on young massive clusters (YMC), that display ideal characteristics to analyse the effect. The majority of them originates from massive star-forming complexes (Grosbøl & Dottori 2012; Kruijssen et al. 2019), where the high densities and the structures created by supersonic motions are likely to have a more prominent impact on their evolution. In addition, their short evolutionary timescale allows us to study the presence of coupling between the processes.

2. Simulation methods

To simulate the interaction between a star cluster and its surrounding environment, each component requires a different approach (N-body and hydrodynamical, respectively). We used specialised codes to follow their evolution and handled the interaction between them with AMUSE (Astrophysical Multipurpose Software Environment; Portegies Zwart et al. 2009; Pelupessy et al. 2013; Portegies Zwart & McMillan 2018). We adopted the N-body code PeTar (Wang et al. 2020), which combines algorithmic regularisations, particle–particle and particle–tree techniques (Kustaanheimo & Stiefel 1965; Ahmad & Cohen 1973; Barnes & Hut 1986), for the cluster dynamics, and the smoothed-particle hydrodynamics (SPH) codes GADGET-2 (Springel 2005) and Fi (Gerritsen & Icke 1997; Pelupessy et al. 2004) to treat the environment. We performed four different sets of simulations: clusters with equal mass stars, clusters with an initial mass function (IMF), and each of these with and without the environment. For convenience, we marked pure N-body reference runs (i.e, without environment) with the suffix Nb. We describe the choice of the initial conditions below.

2.1. Star cluster

We chose the physical parameters to resemble a typical YMC (see, e.g., Portegies Zwart et al. 2010), with the total mass of each model M = 104 M⊙. To study the effects of the interplay of the dynamical processes in the cluster and the environment, we chose five different IMFs: Salpeter (1955) with 0.15 ≤ m/M⊙ ≤ 3 and with 0.15 ≤ m/M⊙ ≤ 6; Kroupa (2001) with 0.1 ≤ m/M⊙ ≤ 3 and with 0.1 ≤ m/M⊙ ≤ 6 (the change from 0.15 to 0.1 M⊙ yields a more similar mean mass to the models with the Salpeter IMF); and an equal-mass stellar population with m = 0.38 M⊙, which is the mean mass of the narrower Salpeter distribution. For the star cluster sizes, we adopted three different values of the virial radius, rvir = 0.7, 1.3 and 3.0 pc. We generated the stellar positions and velocities according to the Plummer (1911) profile. The corresponding core density, ρc, and all combinations of the initial conditions used are listed in Table A.1. The models names are a combination of rvir and the IMF with the upper mass limit (the name contains the suffix “_Nb” when the model evolves only as an isolated N-body system).

The observed relative motions between clusters in star-forming regions are of the order of a few km s−1, leading to an evolution inside the dense region of a few tens of Myr (see Sect. 2.2; Elson et al. 1987; Kuhn et al. 2019; Roccatagliata et al. 2018, 2020). Therefore, we followed the evolution of each star cluster for ≳100 Myr, which would be a conservative choice to let the system escape the high-density environment. We did not include stellar evolution since even the most massive stars in our models (i.e., 6 M⊙) would have the main sequence phase of a similar length as the time span of our simulations (cf. Kippenhahn et al. 2012).

2.2. Environment

We modelled the region surrounding the star cluster with GADGET-2 using equal-mass SPH particles, initially distributed uniformly and at rest in a box of size L = 400 pc, with periodic boundary conditions. The mean number density and the mean molecular weight were the same in all models, that is, ⟨n⟩=10 cm−3 and μ = 1.27, respectively. These values are typical for giant molecular clouds complexes (see, e.g., Blitz & Shu 1980; Elmegreen 2007; Heiner et al. 2008a,b; Ballesteros-Paredes et al. 2020). With the resolution adopted in the reference simulation (1283 SPH particles), these parameters lead to a mass for each gas particles of 9.3 M⊙. In addition, we introduced the gas thermal behaviour according to the method described in Vázquez-Semadeni et al. (2007), adopting the parametric heating and cooling functions from Koyama & Inutsuka (2002).

We gradually developed a divergence-free turbulent field by injecting a fixed amount of energy at each integration time-step, giving a random velocity kick to each gas particle. We defined this energy input using a prescription from Mac Low (1999) such that the equilibrium with the turbulent dissipation would be reached when the box had a velocity dispersion of Vrms (see Table A.2 for the exact values in the simulations). We generate the kicks according to a Gaussian random field having flat spectrum only in between kkick ≤ k ≤ kkick + 1, where we randomly chose the value of kkick at each time-step from a uniform distribution in the range listed in Table A.2.

We did not include the self-gravity of the gas. However, the supersonic turbulent motions (with Vrms ≈ 15 times the sonic speed at equilibrium for gas densities of n = 10 cm−3) and the thermal instability (Field 1965) allowed by the heating and cooling functions are sufficient to create dense structures in the gas. Moreover, in our set-up, we do not expect gravity to play a crucial role on the large scale structures (see also, e.g., Klessen et al. 2000 for an in-depth study of the relative importance of gravity and supersonic motions in highly turbulent environments). Further details of the environmental setup are in Suin (2022).

2.3. Cluster and environment

In the environmental simulations of clusters, we first let the environment fully develop the turbulent field. Once the velocity dispersion of the box stabilises on the desired Vrms, we insert the star cluster. To avoid the phase space of the simulations becoming too wide, we limited this part of the study to a single set of parameters of the background (Vrms = 15 km s−1 and kkick ∈ [4, 8], using a resolution of 1283 SPH particles; see Table A.2). This setup mimics the mean velocity dispersion (Mac Low & Klessen 2004; Dib et al. 2006; Klessen & Hennebelle 2010; Kritsuk et al. 2017) and driving wavelength of the physical regions we are schematising (cf. Norman & Ferrara 1996; Scalo & Elmegreen 2004; Brunt et al. 2009). We also check that there are no massive structures nearby as they could alter the dynamical evolution significantly – we required that the mass of gas included within two virial radii was ≲5% of the cluster mass. The systems then continue to evolve together, with the cluster perceiving the potential generated by the environment. We treated this interaction using the Bridge routine implemented in AMUSE, and letting the code Fi (Gerritsen & Icke 1997; Pelupessy et al. 2004) generate the potential field of the gas (Rieder et al. 2022). We note that the domain of the gaseous environment is periodic. Thus, stars escaping the cluster at higher velocities can depart arbitrarily far from the cluster and still feel the generated potential.

3. Results

We describe the impact of the surrounding environment on the star cluster evolution in three separate ways. First, we focus on the effect of environmental harassment on the core evolution. In the second part, we analyse the cluster outer density distribution. Finally, we look at the overall stellar mass loss.

3.1. Core evolution

Although the tidal force acting on the stars increases linearly with the distance from the cluster centre, the environment also has a profound effect on the inner parts of the cluster – its presence leads to a quantitative acceleration of the core contraction. Tidal shocks increase the velocity dispersion of the cluster so that more stars evaporate from the core to the outer halo. This acts in the same direction as the dynamical relaxation, and the evolution of the core speeds up (Spitzer 1987). We can see this for selected models in Fig. 1 which shows the evolution of the core radius, defined as

r c = ∑ i ρ i 2 r i 2 ∑ i ρ i 2 $$ \begin{aligned} r_{\rm c} = \sqrt{ \frac{ \sum _{i}{\rho ^2_i r^2_i} }{ \sum _{i}{\rho ^2_i} } } \end{aligned} $$(1)

(cf. Aarseth 2003). Here ρi is the density estimator near the ith cluster member, computed with the Hop clump-find algorithm (Eisenstein & Hut 1998) using its nearest 12 neighbours, and ri is its distance from the cluster density centre (Casertano & Hut 1985). Both compact models with an IMF experience core collapse (highlighted with a vertical line in the left-hand plots of Fig. 1), but the environmental run reaches it sooner. This faster contraction of the core is also evident in the models where core collapse did not happen within the simulated temporal interval (see how the lines spread in the right-hand plot). Moreover, the plots show that the presence of a background has almost no effect on the core evolution of the equal mass clusters (the black lines remains superimposed during the whole simulation). The same is true during the very first phases of evolution in cluster with IMF (up to 10 Myr in the compact cluster and 30 Myr in the other).

thumbnail Fig. 1.

Evolution of the core radius, see Eq. (1), normalised to the initial core radius. The plots show the clusters with the Salpeter IMF or equal-mass stars, both isolated and embedded into the environment (as labelled in the legends, we note that the two black lines are almost superimposed). The left-hand plot shows the more compact models with rvir = 0.7 pc, the right-hand one is for rvir = 1.3 pc. The moments of core collapse are highlighted with a vertical line.

Although useful for tracking the temporal development of the core, Eq. (1) is a functional form that differs from the actual core radius in a way that depends on the shape of the cluster itself (see, e.g., Table 2 in Casertano & Hut 1985). To identify the moment of core collapse, tcc, we followed the self-similar argument of Lynden-Bell & Eggleton (1980). The core radius evolves as

r c ( t ) ∝ ( t cc − t ) 2 / ( 6 − α ) , $$ \begin{aligned} r_{\rm c}(t) \propto \left( t_\mathrm{cc} -t \right)^{2/(6-\alpha )}\!, \end{aligned} $$(2)

for t ≤ tcc, where α is a fitting parameter determined through numerical simulations. To apply this equation, we used the procedure of Pavlík & Šubr (2018) who show that before core collapse, the minima of Lagrangian radii are a fixed multiple of rc. Consequently, the evolution of these minima follows Eq. (2)1.

Table 1 displays the results of the fit for the core collapsed clusters. In all simulations, we recover the expected value of α ≈ 2.21 (e.g., Takahashi 1995; Lynden-Bell & Eggleton 1980) within 3 σ uncertainties. However, there is a tendency towards lower values in the simulations with an upper mass limit of 6 M⊙. This is compatible with the result of Pavlík & Šubr (2018), who found α ≈ 1.5 using a Salpeter IMF with an about ten times wider range than ours. The R13_Sal6_Nb model shows the greatest error in both fit parameters. In this simulation, the core contraction is less pronounced, which makes the location of the Lagrangian radii minima extremely sensitive to random fluctuations. Concerning the moment of collapse, Table 1 shows that the gaseous background accelerates the cluster evolution. Environmental simulations display smaller tcc with respect to their N-body counterpart. The only exception is run R07_Sal6, for which the dynamical timescale is short and the effect of the environment is negligible.

Table 1.

List of fit results.

We also analysed the binary ejected by the clusters after core collapse. There is a good correlation between the formation of the first tight binary (cf. Fujii & Zwart 2014; Pavlík & Šubr 2018). We repeated the simulation R07_Sal6 four times with the environment and eight without increasing the size of the ejected binaries sample. We did not detect any statistically significant difference between the binding energy distributions extracted from the two cases. We refer to Suin (2022) for the in-depth study of the binary characteristics.

3.2. Outer region

We evaluated the impact of the background on the evolution of the radial density profile of the outer region. We chose the radial shells between the 94% and 98% Lagrangian radii with binning by 0.5% and performed a least-square fit for the power-law

ρ ( r ) ∝ r − β . $$ \begin{aligned} \rho (r)\propto r^{-\beta }. \end{aligned} $$(3)

We selected the lower bound to achieve an error smaller than 5% for a Plummer model fit2. The upper bound reduces the impact of statistical relative fluctuations (bigger in the more sparse regions) on the fit. The procedure implicitly assumes spherical symmetry. This, however, is justified by Fig. 2, which shows the behaviour of the surface density in run R07_Sal3 projected on the planes xz and yz. The two sequences closely follow each other, so the three-dimensional nature of the problem only marginally affects the exponent β.

thumbnail Fig. 2.

Temporal evolution of the modulus of the surface density exponent in the simulation R07_Sal3, projected onto the planes xz (blue) and yz (red). We smoothed the data with a simple moving weighted average between 5 points (2.5 Myr). The black line shows the value of Ξtid (axes on the right), computed at the centre of the cluster.

To compare the behaviour of β with the tidal field acting on the cluster, we defined the variable Ξtid(t) as the sum in quadrature of tidal tensor elements, computed at the cluster core centre xc,

Ξ tid ( t ) ≡ Ξ α β ( x c , t ) Ξ α β ( x c , t ) , $$ \begin{aligned} \Xi _\mathrm{tid} (t) \equiv \sqrt{\Xi _{\alpha \beta }(\boldsymbol{x_\mathrm{c} }, t) \Xi _{\alpha \beta }(\boldsymbol{x_\mathrm{c} }, t)}, \end{aligned} $$(4)

where

Ξ α β ( x c , t ) ≡ − ∂ 2 ϕ ∂ x α ∂ x β ( x c , t ) , $$ \begin{aligned} \Xi _{\alpha \beta } (\boldsymbol{x_\mathrm{c} }, t) \equiv -\frac{\partial ^2\phi }{\partial x_\alpha \partial x_\beta }(\boldsymbol{x_\mathrm{c} }, t), \end{aligned} $$(5)

and ϕ(x, t) is the gravitational potential per unit mass generated by the gas at time t. We adopt the Einstein convention, summing over repeated indices. Ξtid(t) has dimension of Myr−2 and is invariant under the change of coordinates, so it provides a good indicator of the strength of the tidal field affecting the cluster. To give an idea of the physical value of Ξtid, a point-like object of mass Mg and distance Rg from the cluster would exert Ξ tid = 6 G M g / R g 3 $ \Xi_{\mathrm{tid}}=\sqrt{6}\, G M_{\mathrm{g}}/R_{\mathrm{g}}^{3} $. In physical units, this scales as

Ξ tid ∼ 0.01 ( M g 10 3 M ⊙ ) ( R g 10 pc ) − 3 Myr − 2 . $$ \begin{aligned} \Xi _\mathrm{tid} \sim 0.01\left(\frac{M_\mathrm{g} }{10^3\, M_\odot }\right)\left(\frac{R_\mathrm{g} }{10\, \mathrm{pc}}\right)^{-3}\, \mathrm{Myr}^{-2}. \end{aligned} $$(6)

The trace, in this case, would be zero.

Figure 3 shows the temporal evolution of β in the simulations with and without the ambient gas, along with the indicator Ξtid of the environmental run. In the isolated simulations, β remains almost constant initially, fluctuating around the initial value of 5 (which comes from the Plummer model). It then starts to decrease approximately linearly in time while approaching core collapse.

thumbnail Fig. 3.

Evolution of β throughout a selected set of simulations, with the axis on the left. The shaded areas identify the value of the fit within 1σ (blue for the N-body runs, red for those with environment). We smoothed the data with a simple moving weighted average of width 2.5 Myr (5 model outputs). The black line shows the value of Ξtid of the environmental runs, as in Fig. 2 (axis on the right). The vertical dashed lines in the plots (b) and (e) highlight the moment of core collapse (see also Table 1). (a) R07_NoImf, (b) R07_Sal3, (c) R13_NoImf, (d) R13_Sal3, (e) R13_Sal6, (f) R30_Sal6.

Runs with the environment show greater fluctuations of β. In addition, we see that the effect of the background increases if the cluster is wider and its evolutionary phase more advanced. Large clusters are more sensitive to tides, and the core contraction sends stars out to the halo, where tides are stronger. A comparison between R07_NoImf (Fig. 3a) and R13_NoImf (Fig. 3c) allows us to discuss the first aspect. Without a mass function, the dynamical evolution of these clusters is almost negligible throughout the simulations (see the black lines in Fig. 1), and we can extract the effect of the tides on the evolution directly. The clusters do not significantly differ from the evolution of their isolated counterpart. However, the fluctuations are bigger in both cases. The same can also be inferred from R30_Sal6 (Fig. 3f) – in this case, the low stellar density slows down the mass segregation process and makes its dynamical development similarly slow as in the equal-mass case. In addition, in R13_NoImf (Fig. 3c) the trigger for the largest fluctuations in β is the major encounter at ≈10,  40 , 65 Myr, although in run R07_NoImf it is not possible to extract a clear correspondence between close encounters (peaks in Ξtid) and variation of the power-law slope.

Widening the mass range decreases the evolutionary timescale, and enhances the effect of the background. Core collapsed clusters (Figs. 3a–e) display a more step-like evolution of the exponent β when inserted into the environment. While approaching core collapse, the density gradient becomes more sensitive to the peaks in Ξtid. Close to these peaks, β initially increases and the cluster becomes more compact. Then, stars populate the outer halo and β decreases. The situation is similar in Gieles et al. (2006), where the simplest example consists in a head-on encounter of the cluster with a spherical cloud. While approaching the cloud, the cluster contracts as stars suffer a net acceleration toward the centre along the direction perpendicular to the motion. After the encounter, the cluster expands due to the relative perpendicular velocities acquired by particles. The parallel velocity variation, instead, cancels out during the whole encounter.

By the end of the simulation, the harassed clusters display softer outer power-law slopes. This also shows up in the last ≈20 Myr of model R13_Sal3 (Fig. 3d), which at the end of the simulation is close to experiencing core collapse (as can also be deduced from Fig. 1). When far from core collapse, these clusters behave similarly to the equal mass case. At early stages of evolution, its shape did not change, and the environment does not affect the outer region permanently. Even though showing larger fluctuation than the reference runs, β does not depart significantly from the initial value. As soon as a tidal shock occurs and the cluster has an evolved configuration, β departs from the linear behaviour of the isolated case (see the shocks at ≈40,  70 Myr and the subsequent decrease in β in Figs. 3a–e).

3.3. Cluster mass loss

To determine whether a star is bound to the cluster is not straightforward when an external potential is present. The cluster members are constantly interacting with the gaseous environment which can differentially accelerate them. While passing through or near a cloud, a star in the cluster gains enough kinetic energy to appear unbound, and the usual definition, which labels escaping stars based on their instantaneous binding energy to the system, fails. However, over a whole encounter, only the orthogonal velocity component remains high, while any velocity gains in the parallel component cancel out. To avoid the fluctuations in the mass loss due to these encounters, we only mark as “escapers” the stars that remain unbound throughout the last 5 Myr of the simulation (see Appendix B for further justification of this choice).

Figure 4 shows two examples of the cumulative mass loss over time, comparing the simulations with and without the environment, along with the mass loss from the equal mass clusters. This highlights the presence of a coupling mechanism between the two evolutions. As for the outer density distribution, the coupling only shows up when both effects are present. The mass loss of run R07_Sal3 is particularly explanatory (see the top panel of Fig. 4), as it follows the isolated cluster until just before core collapse. In addition, after core collapse the mass loss presents a steeper increase than in the isolated cluster.

thumbnail Fig. 4.

Evolution of the cluster mass loss in percents of the total initial mass (see the text for the definition).

Table 2 provides the mass lost in each simulation after 95 Myr. In every simulation, the mass loss from the isolated equal mass cluster is negligible. This supports the choice of using the embedded equal mass cluster as a reference for the evolution caused only by environmental harassment. As expected, the importance of the tides increases as the cluster expands. In contrast, the mass lost in pure N-body simulations depends only on the mass segregation timescale, which is proportional to ρ−1/2 and inversely proportional to the width of the IMF.

Table 2.

Mass loss of simulated cluster computed at 95 Myr.

Finally, it is worth noting that two-body relaxation and tidal interactions affect the IMF in different ways. The former causes mass segregation with higher-mass stars moving inward and lower-mass stars outward, which enables the lower-mass stars to evaporate from the cluster (see also, e.g., Chandrasekhar 1942; Spitzer 1969; Binney & Tremaine 2008). The cluster then forms a dense core of high-mass stars and once it starts to collapse, binaries form in the centre. These perturb their neighbours, and we begin to observe massive stars escaping from the mass segregated core (see also Hills 1975; Hut 1985). On the other hand, if the relaxation timescales are long, the cluster does not segregate as quickly and the tidal harassment is the dominant reason for mass-loss. Stars are peeled-off from the outer regions independent of their masses. Consequently, the relative contributions of these two processes determine the mass distribution of the escapers. This is shown in Fig. 5, where the mass distribution of escaping stars divided by the original IMF in models R07_Sal13 and R07_Sal13_Nb. When the environment is present, the tidal shocks enhance fraction of escaping stars at intermediate masses. Instead, the loss of massive stars, similar in the two models, can be linked to the relaxation process ongoing in the core that is composed almost exclusively of massive stars and binaries at later times. In our model, the percentage of mass loss is small so this does not significantly affect the mass function of the remaining members. However, the same process happens in protoclusters still embedded in the parent cloud. Their gas density is orders of magnitude higher than our n = 10 cm−3, and structures are closer to the cluster (Kruijssen et al. 2012; Kruijssen 2012). At the same time, the IMF of massive protoclusters is much wider, since massive stars are still on the main sequence, so that the evolutionary timescale is also reduced (Allison et al. 2009; Yu et al. 2011). In the end, the effects we described could indeed play a crucial role when inferring the IMF of young clusters.

thumbnail Fig. 5.

Mass distribution of escapers normalised to the initial IMF for the runs R07_Sal3 (red) and its N-body reference run (blue).

4. Limitations

The choice of initial conditions has trade-offs. Adopting the Plummer model and truncating the IMF at 6 M⊙ allows us to avoid the more complex stages of the cluster formation that a more physical model would require. However, stellar evolution and feedback play a significant role when including younger stars with higher masses. Massive stars sink faster toward the cluster centre and die earlier (Portegies Zwart et al. 2007). The presence of a denser environment also enhances tidal perturbations. Nevertheless, the mechanisms studied here will impact the evolution of the cluster. Regarding the background setup, we mimicked the averaged observed values of star-forming regions. The constant energy input generated a constant turbulent field in which structures formed and dissolved. The feedback mechanisms acting in such environments are complicated and still far from completely understood. Adding self-gravity increases the computational times due to the collapsing high-density regions. Our work is, instead, intended to show that even simple models have observational consequences and to provide insights regarding the interplay between the dynamical evolution of star clusters and the gaseous background that will only be magnified by achieving more physical simulations.

5. Conclusions

We conducted a series of numerical simulations following the evolution of both the star cluster and its surrounding ISM environment. Gaseous structures, created by supersonically driven turbulent motions, produced a spatially and temporally varying external potential which perturbed the embedded cluster.

The core evolution and the external power-law slope of the density distribution both exhibited effects induced by the gas. The tidal shocks increase the fractions of core stars that are sent into the halo. This process acts similarly to dynamical relaxation. We note that this sped-up redistribution of energy can also continue to later stages of the cluster evolution, every time the cluster passes through a star-forming region. This supports the finding of Gnedin et al. (1999) of a reduced core collapse timescale using spherical Fokker–Planck models.

For the cluster periphery, we found that the exponent of the asymptotic density distribution, β, displays greater fluctuations induced by the environment. While approaching core collapse, tides significantly affect β. The slope departs from the almost linear temporal evolution expected in the isolated case, displaying a more step-like behaviour. In addition, all core collapsed runs end up with a lower value of β as a result of the accelerated evolution. Both effects would contribute to observing more clusters in a more relaxed and expanded configuration (e.g., Portegies Zwart et al. 2010; Banerjee & Kroupa 2017).

This coupling also manifests in the cluster mass loss. In each environmental simulation, the final mass loss is larger than the sum of the respective reference runs. This highlights a coupling between the dynamical evolution and the presence of the environment.

1

During their evolution, Lagrangian radii exhibit large fluctuations. These are particularly important in the inner regions, where the number of star is small, leading to a high statistical noise. To achieve a greater accuracy in the localisation of the minima, we smoothed the Lagrangian radii, using Savitzky & Golay (1964) smoothing algorithm with a second degree polynomial, as in Pavlík & Šubr (2018). The smoothing window corresponded to 9 data-points, i.e., 4.5 Myr. The algorithm is implemented in Python as scipy.signal.savgol_filter.

2

For a Plummer model of scale radius a, β = 5x2/(1 + x2), where x ≡ r/a. Taking r to be the 94% Lagrangian radius, the exponent differs from its asymptotic value by ≈0.045.

Acknowledgments

We thank Maurizio Davini and the Direzione Infrastrutture Digitali of the University of Pisa for providing access to two 128-core DELL R630i servers of the High-Performance Computing Division of the San Piero a Grado Green Data Center, without which this work would not have been possible. We are grateful to Veronica Roccatagliata and Michele Cignoni for insightful discussions. We also thank the AMUSE team for the support provided. This research has made use of NASA’s Astrophysics Data System Bibliographic Services. The Python programming language with NumPy (Harris et al. 2020) and Matplotlib (Hunter 2007) were used in this project.

References

  1. Aarseth, S. J. 2003, Cambridge: Gravitational N-Body Simulations (Cambridge: University of Cambridge Press) [CrossRef] [Google Scholar]
  2. Adamo, A., Zeidler, P., Kruijssen, J. M., et al. 2020, Space Sci. Rev., 216, 69 [NASA ADS] [CrossRef] [Google Scholar]
  3. Ahmad, A., & Cohen, L. 1973, J. Comput. Phys., 12, 389 [Google Scholar]
  4. Allison, R. J., Goodwin, S. P., Parker, R. J., et al. 2009, ApJ, 700, L99 [Google Scholar]
  5. Ballesteros-Paredes, J., André, P., Hennebelle, P., et al. 2020, Space Sci. Rev., 216, 76 [Google Scholar]
  6. Banerjee, S., & Kroupa, P. 2017, A&A, 597, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Barnes, J., & Hut, P. 1986, Nature, 324, 446 [NASA ADS] [CrossRef] [Google Scholar]
  8. Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton: Princeton University Press) [Google Scholar]
  9. Blitz, L., & Shu, F. H. 1980, ApJ, 238, 148 [NASA ADS] [CrossRef] [Google Scholar]
  10. Brunt, C. M., Heyer, M. H., & Mac Low, M. M. 2009, A&A, 504, 883 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Casertano, S., & Hut, P. 1985, ApJ, 298, 80 [Google Scholar]
  12. Chandrasekhar, S. 1942, Principles of Stellar Dynamics (Chicago: University of Chicago Press), 321 [Google Scholar]
  13. Chernoff, D. F., & Weinberg, M. D. 1990, ApJ, 351, 121 [NASA ADS] [CrossRef] [Google Scholar]
  14. Chevance, M., Kruijssen, J. M. D., Vazquez-Semadeni, E., et al. 2020, Space Sci. Rev., 216, 50 [CrossRef] [Google Scholar]
  15. Colín, P., Vázquez-Semadeni, E., & Gómez, G. C. 2013, MNRAS, 435, 1701 [CrossRef] [Google Scholar]
  16. Dib, S., Bell, E., & Burkert, A. 2006, ApJ, 638, 797 [NASA ADS] [CrossRef] [Google Scholar]
  17. Eisenstein, D. J., & Hut, P. 1998, ApJ, 498, 137 [NASA ADS] [CrossRef] [Google Scholar]
  18. Elmegreen, B. G. 2007, ApJ, 668, 1064 [NASA ADS] [CrossRef] [Google Scholar]
  19. Elson, R. A. W., Fall, S. M., & Freeman, K. C. 1987, ApJ, 323, 54 [Google Scholar]
  20. Federrath, C., & Klessen, R. S. 2013, ApJ, 763, 51 [Google Scholar]
  21. Field, G. B. 1965, ApJ, 142, 531 [Google Scholar]
  22. Fujii, M. S., & Zwart, S. P. 2014, MNRAS, 439, 1003 [NASA ADS] [CrossRef] [Google Scholar]
  23. Gerritsen, J. P., & Icke, V. 1997, A&A, 325, 972 [NASA ADS] [Google Scholar]
  24. Gieles, M., Portegies Zwart, S. F., Baumgardt, H., et al. 2006, MNRAS, 371, 793 [Google Scholar]
  25. Gnedin, O. Y., & Ostriker, J. P. 1997, ApJ, 474, 223 [Google Scholar]
  26. Gnedin, O. Y., Lee, H. M., & Ostriker, J. P. 1999, ApJ, 522, 935 [Google Scholar]
  27. Grosbøl, P., & Dottori, H. 2012, A&A, 542, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357 [Google Scholar]
  29. Heiner, J. S., Allen, R. J., Emonts, B. H. C., & van der Kruit, P. C. 2008a, ApJ, 673, 798 [NASA ADS] [CrossRef] [Google Scholar]
  30. Heiner, J. S., Allen, R. J., Wong, O. I., & Van Der Kruit, P. C. 2008b, A&A, 489, 533 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  31. Hills, J. G. 1975, AJ, 80, 809 [NASA ADS] [CrossRef] [Google Scholar]
  32. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  33. Hut, P. 1985, in Dynamics of Star Clusters: Proceedings of the 113th Symposium of the International Astronomical Union, eds. J. Goodman, & P. Hut (Princeton: Reidel Publishing Company), 231 [Google Scholar]
  34. Kim, J., Chevance, M., Diederik Kruijssen, J. M., et al. 2021, MNRAS, 504, 487 [NASA ADS] [CrossRef] [Google Scholar]
  35. Kippenhahn, R., Weigert, A., & Weiss, A. 2012, Stellar Structure and Evolution, Astronomy and Astrophysics Library (Berlin: Springer) [CrossRef] [Google Scholar]
  36. Klessen, R. S., & Hennebelle, P. 2010, A&A, 520, A17 [NASA ADS] [CrossRef] [Google Scholar]
  37. Klessen, R. S., Heitsch, F., & Mac Low, M.-M. 2000, ApJ, 535, 887 [NASA ADS] [CrossRef] [Google Scholar]
  38. Koyama, H., & Inutsuka, S.-I. 2002, ApJ, 564, L97 [Google Scholar]
  39. Kritsuk, A. G., Ustyugov, S. D., & Norman, M. L. 2017, New J. Phys., 19, 065003 [Google Scholar]
  40. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  41. Kruijssen, J. M. D. 2012, MNRAS, 426, 3008 [Google Scholar]
  42. Kruijssen, J. M. D., Maschberger, T., Moeckel, N., et al. 2012, MNRAS, 419, 841 [NASA ADS] [CrossRef] [Google Scholar]
  43. Kruijssen, J. M. D., Schruba, A., Chevance, M., et al. 2019, Nature, 569, 519 [NASA ADS] [CrossRef] [Google Scholar]
  44. Kuhn, M. A., Hillenbrand, L. A., Sills, A., Feigelson, E. D., & Getman, K. V. 2019, ApJ, 870, 32 [CrossRef] [Google Scholar]
  45. Kustaanheimo, P., & Stiefel, E. 1965, J. Reine Angew. Math., 1965, 204 [CrossRef] [Google Scholar]
  46. Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57 [Google Scholar]
  47. Lynden-Bell, D., & Eggleton, P. P. 1980, MNRAS, 191, 483 [NASA ADS] [Google Scholar]
  48. Mac Low, M. 1999, ApJ, 524, 169 [NASA ADS] [CrossRef] [Google Scholar]
  49. Mac Low, M.-M., & Klessen, R. S. 2004, Rev. Mod. Phys., 76, 125 [Google Scholar]
  50. Norman, C. A., & Ferrara, A. 1996, ApJ, 467, 280 [NASA ADS] [CrossRef] [Google Scholar]
  51. Pavlík, V., & Šubr, L. 2018, A&A, 620, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  52. Pelupessy, F. I., Van Der Werf, P. P., & Icke, V. 2004, A&A, 422, 55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Pelupessy, F. I., van Elteren, A., de Vries, N., et al. 2013, A&A, 557, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Plummer, H. C. 1911, MNRAS, 71, 460 [Google Scholar]
  55. Portegies Zwart, S., & McMillan, S. 2018, Astrophysical Recipes: The art of AMUSE (IOP Publishing) [CrossRef] [Google Scholar]
  56. Portegies Zwart, S. F., McMillan, S. L., & Makino, J. 2007, MNRAS, 374, 95 [NASA ADS] [CrossRef] [Google Scholar]
  57. Portegies Zwart, S., McMillan, S., Harfst, S., et al. 2009, New Astron., 14, 369 [Google Scholar]
  58. Portegies Zwart, S. F., McMillan, S. L., & Gieles, M. 2010, ARA&A, 48, 431 [NASA ADS] [CrossRef] [Google Scholar]
  59. Rieder, S., Dobbs, C., Bending, T., Liow, K. Y., & Wurster, J. 2022, MNRAS, 509, 6155 [Google Scholar]
  60. Roccatagliata, V., Sacco, G. G., Franciosini, E., & Randich, S. 2018, A&A, 617, L4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  61. Roccatagliata, V., Franciosini, E., Sacco, G. G., Randich, S., & Sicilia-Aguilar, A. 2020, A&A, 638, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  62. Salpeter, E. E. 1955, ApJ, 121, 161 [Google Scholar]
  63. Savitzky, A., & Golay, M. 1964, J. Anal. Chem., 36, 1627 [Google Scholar]
  64. Scalo, J., & Elmegreen, B. G. 2004, ARA&A, 42, 275 [NASA ADS] [CrossRef] [Google Scholar]
  65. Spitzer, L. 1969, ApJ, 158, L139 [NASA ADS] [CrossRef] [Google Scholar]
  66. Spitzer, L. 1987, Dynamical Evolution of Globular Clusters (Princeton: University of Princeton Press) [Google Scholar]
  67. Spitzer, L., & Shapley, H. 1940, MNRAS, 100, 396 [NASA ADS] [CrossRef] [Google Scholar]
  68. Springel, V. 2005, MNRAS, 364, 1105 [Google Scholar]
  69. Suin, P. 2022, Master’s Thesis, Department of Physics, University of Pisa, Italy [Google Scholar]
  70. Takahashi, K. 1995, PASJ, 47, 561 [Google Scholar]
  71. Terlevich, E. 1987, MNRAS, 224, 193 [NASA ADS] [CrossRef] [Google Scholar]
  72. Vázquez-Semadeni, E. 1994, ApJ, 423, 681 [CrossRef] [Google Scholar]
  73. Vázquez-Semadeni, E., Gómez, G. C., Jappsen, A. K., et al. 2007, ApJ, 657, 870 [CrossRef] [Google Scholar]
  74. Wang, L., Iwasawa, M., Nitadori, K., & Makino, J. 2020, MNRAS, 497, 536 [NASA ADS] [CrossRef] [Google Scholar]
  75. Yu, J., de Grijs, R., & Chen, L. 2011, ApJ, 732, 16 [Google Scholar]

Appendix A: Model parameters

Table A.1.

List of N-body simulations of the star cluster.

Table A.2.

List of SPH simulations of the environment.

Appendix B: Stellar mass loss

Encounters with the environment increase the velocities of stars relative to the cluster centre, making them appear unbound for a short period of time (see the red line in Fig. B.1. The spikes, which show the instantaneous increase of the number of stars with positive energy, correspond to the moments when the cluster passes near a cloud. After the passage, most stars decelerate and return to a bound state. To minimise these episodic fluctuations, we consider as “escapers” only stars that appear unbound for the last 5 Myr of the simulation (see the black lines in Fig. B.1, and also Mlost in Sect. 3.2, Fig. 4 and Tab. 2). The choice of this particular interval is dictated by the time with which the autocorrelation of Ξtid, averaged between the simulations, crosses zero (see Suin 2022 for further details). An alternative method is to label as escapers those stars that remain unbound for at least some predefined interval, τw, starting from the time of their first detection as unbound. In this way, the evolution of Mlost is independent of the final state of the system. A limitation is that its monotonicity might be lost if τw were too small (see, e.g., the yellow line in Fig. B.1). If, however, τw = 5 Myr this method and the one we used in Sect. 3.2 are equivalent, with the latter being less computationally demanding.

thumbnail Fig. B.1.

Comparison of the methods to compute the fractional mass loss in model R07_Sal6. Left: Comparison between our model (black) and the mass loss calculated instantaneously (red). Right: The alternative method proposed in the text compared to our model (black line), coloured lines correspond to different window lengths τw.

All Tables

Table 1.

List of fit results.

In the text
Table 2.

Mass loss of simulated cluster computed at 95 Myr.

In the text
Table A.1.

List of N-body simulations of the star cluster.

In the text
Table A.2.

List of SPH simulations of the environment.

In the text

All Figures

thumbnail Fig. 1.

Evolution of the core radius, see Eq. (1), normalised to the initial core radius. The plots show the clusters with the Salpeter IMF or equal-mass stars, both isolated and embedded into the environment (as labelled in the legends, we note that the two black lines are almost superimposed). The left-hand plot shows the more compact models with rvir = 0.7 pc, the right-hand one is for rvir = 1.3 pc. The moments of core collapse are highlighted with a vertical line.
In the text
thumbnail Fig. 2.

Temporal evolution of the modulus of the surface density exponent in the simulation R07_Sal3, projected onto the planes xz (blue) and yz (red). We smoothed the data with a simple moving weighted average between 5 points (2.5 Myr). The black line shows the value of Ξtid (axes on the right), computed at the centre of the cluster.
In the text
thumbnail Fig. 3.

Evolution of β throughout a selected set of simulations, with the axis on the left. The shaded areas identify the value of the fit within 1σ (blue for the N-body runs, red for those with environment). We smoothed the data with a simple moving weighted average of width 2.5 Myr (5 model outputs). The black line shows the value of Ξtid of the environmental runs, as in Fig. 2 (axis on the right). The vertical dashed lines in the plots (b) and (e) highlight the moment of core collapse (see also Table 1). (a) R07_NoImf, (b) R07_Sal3, (c) R13_NoImf, (d) R13_Sal3, (e) R13_Sal6, (f) R30_Sal6.
In the text
thumbnail Fig. 4.

Evolution of the cluster mass loss in percents of the total initial mass (see the text for the definition).
In the text
thumbnail Fig. 5.

Mass distribution of escapers normalised to the initial IMF for the runs R07_Sal3 (red) and its N-body reference run (blue).
In the text
thumbnail Fig. B.1.

Comparison of the methods to compute the fractional mass loss in model R07_Sal6. Left: Comparison between our model (black) and the mass loss calculated instantaneously (red). Right: The alternative method proposed in the text compared to our model (black line), coloured lines correspond to different window lengths τw.
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.