© ESO 2020

1. Introduction

The discovery of supermassive black holes (SMBHs) in the first billion years of the Universe (Fan et al. 2006; Mortlock et al. 2011; Wu et al. 2015; Reed et al. 2019) has led to the study of formation channels for massive black hole seeds early on (z ∼ 10−20), following the collapse of molecular-cooling halos (∼10⁶ M_⊙) or atomic-cooling halos (∼10⁸ M_⊙). The direct collapse of protogalactic clouds exposed to a moderate Lyman-Werner flux from neighboring star forming halos produced the most massive seeds (10³⁻⁵ M_⊙) (Wise et al. 2019), and recent numerical simulations suggest that contrary to what was previously thought, the UV flux is not that critical anymore in the case of the efficient merger of fragments (Suazo et al. 2019). Another promising mechanism to form massive seeds is the formation of very massive stars (VMS) in the centers of dense stellar systems through stellar mergers (Fujii & Portegies Zwart 2013; Katz et al. 2015; Sakurai et al. 2017; Reinoso et al. 2018), producing either VMSs with several 10² M_⊙ or intermediate mass black holes with ∼10³ M_⊙ and potentially up to 10⁴ M_⊙. Such massive seeds could be present in the early Universe if accreting population III (Pop. III) protostars can merge before entering the main sequence or if the remaining gas in the cluster can be accreted by the central object (Boekholt et al. 2018). Furthermore, after the formation of a BH in the center of a stellar cluster, additional growth can be expected by tidal disruption events of stars passing close to the black hole (Sakurai et al. 2019; Bonetti et al. 2020).

Although most of these previous studies have focused on mergers in second generation star clusters, that is, clusters that formed from molecular clouds that were polluted by stellar winds or supernovae from the first generation of stars (Katz et al. 2015; Sakurai et al. 2017, 2019), we are mostly interested in the very first star clusters of the Universe given the particular properties of these systems and these stars. Fragmentation occurs at higher densities in primordial gas clouds (10⁹ cm⁻³ or higher) (e.g., Clark et al. 2011a; Greif et al. 2011, 2012; Smith et al. 2011, 2012; Latif et al. 2013a), with clusters potentially having half-mass radii of 0.1 pc, especially if dust grains are present in such clouds which trigger fragmentation at high densities (e.g., Omukai et al. 2005; Klessen et al. 2012; Bovino et al. 2016; Latif et al. 2016).

In those environments, the protostellar radii could also be enhanced due to the rapid accretion expected in primordial or low-metallicity gas, which in turn are a consequence of much higher gas temperatures than at present day star formation, and so the protostellar radii are enhanced up to 300 R_⊙ (Stahler et al. 1986; Omukai & Palla 2001, 2003). Moreover, high accretion rates of ∼0.1 M_⊙ yr⁻¹ have been reported in several simulations for these protostars (Hosokawa et al. 2012, 2013; Schleicher et al. 2013; Haemmerlé et al. 2018; Woods et al. 2017), implying ∼500 R_⊙ for a 10 M_⊙ star and potentially more than 1000 R_⊙ for a 100 M_⊙ star. While some models already included the effects of an external potential in the evolution of dense star clusters (e.g., Leigh et al. 2013a,b, 2014; Boekholt et al. 2018; Sakurai et al. 2019) in order to mimic the effects of gas during the embedded phase of the cluster or to account for the dark matter halo, disentangling the effects of the external potential can be tricky because these models also include mass accretion recipes (Boekholt et al. 2018).

In this paper, we present a systematic investigation on how the formation of very massive stars in dense star clusters and the evolution of the clusters themselves are affected when including a background potential in the evolution of the systems. Our work can be used to constraint the parameter space that future, more sophisticated simulations should focus on. Since each physical collision led to a merger in our runs, we use these terms interchangeably. We describe our setup and initial conditions in Sect. 2 and our results are presented in Sect. 3, including the typical evolution of clusters under the influence of a background potential and clusters without a background potential, the number of collisions, the mass of the resulting object, and our model to estimate the number of mergers at different times. We apply our model to population III (Pop.III) star clusters in Sect. 4. A discussion about neglected effects and considerations for future research on this subject is given in Sect. 5 and finally the conclusions are presented in Sect. 6.

2. Simulation setup

To understand the impact of the gas potential on the number of mergers that occur in star clusters, we first performed simulations in which we did not include stellar mergers, and we compared the evolution of these systems without an external potential and when placed in the center of a background potential. Once we addressed the effects of the extra force on the evolution of the clusters, we performed simulations that include stellar mergers in clusters with and without an external potential. To perform the calculations, we used a modified version of NBODY6¹ (Aarseth et al. 2000) to treat collisions, where we switched off the stellar evolution package and instead explicitly specify the stellar radii to perform a parameter study. NBODY6 is a fourth order Hermite integrator, and includes a spatial hierarchy to speed up the calculations: This is referred to as the Ahmad-Cohen scheme (Ahmad & Cohen 1973). It also includes routines to treat tidal capture and tidal circularization of binary systems (Mardling & Aarseth 2001) that we activated in our simulations. Another important routine included is the Kustaanheimo-Stiefel regularization (Kustaanheimo & Stiefel 1965), which is an algorithm that can be used to treat binaries and close two-body encounters more accurately. We performed a total of 344 N-body simulations.

2.1. Simulations without stellar mergers

In order to create a controlled experiment to which we can compare our results including stellar mergers, we first investigated the effects produced by the external potential on the evolution of the star clusters in the absence of collisions. We specifically explored how the mass of the external potential affects the core-collapse, the formation of binary systems, and the evaporation of the clusters. We modeled each cluster as a Plummer (Plummer 1911) distribution of N = 1000 and N = 10 000 equal mass particles with a total cluster mass of M_stars = 10⁴ M_⊙. Given that we are interested in the overall evolution of the system, and in order to save computational time, we modeled less dense clusters than the ones presented in Sect. 2.2, with a virial radius of R_v = 1.0 pc. We then included an analytic background potential that also follows a Plummer density profile with the same virial radius as the stellar distribution R_v, ext = 1.0 pc and the mass of the potential was varied as M_ext = 0.0, 0.1, 0.5 and 1.0 × M_stars. The clusters start in virial equilibrium. We performed a total of eight simulations that do not include stellar mergers and which are listed in Table 1.

2.2. Simulations that include stellar mergers

Our simulations that include stellar mergers aim to model the first stages after the formation of Pop. III star clusters. We modeled a compact cluster in virial equilibrium by also using a Plummer distribution (Plummer 1911) for the stars that are all equal in mass and radius at the beginning of the simulation. We modeled dense clusters with a virial radius of R_v = 0.14 pc and a total mass in stars of M_stars = 10⁴ M_⊙. We then performed the same set of simulations including an external analytic Plummer potential with M_ext = 1.0 × M_stars with the same virial radius of R_v, ext = 0.14 pc in order to consider at first order the effects of the gas that remains in the cluster after the formation of the stars. Taking the mass of the external potential into account, the crossing time of the clusters were calculated as:

$$\begin{array}{c}\hfill t_{\mathrm{cross}}=\sqrt{\frac{R_{\mathrm{v}}^3}{G{M}_{\mathrm{stars}}}}\frac{1}{1+q},\end{array}$$(1)

with $q=\frac{M_{\mathrm{ext}}}{M_{\mathrm{stars}}}$ (see Appendix A); this gives a value of t_cross = 0.0078 Myr and t_cross, ext = 0.0039 Myr for clusters without and with an external potential, respectively. We investigated how the mass of the final merger product depends on the initial number of stars in the cluster N, which we varied as N = 100, 500, 1000, and 5000 stars. We also varied the initial stellar radius as R_star = 20, 50, 100, 200, 500, 1000, and 5000 R_⊙ using equal radii stars in each run. For each of these configurations, we performed simulations with M_ext = 0.0 and 1.0× M_stars. Finally, we ran a total of six random simulations per each configuration, which are listed in Table 2. This gives a total of 336 N-body simulations that include stellar mergers.

2.3. Stellar mergers

In our simulations, a merger occurs when two stars are separated by a distance equal to or smaller than the sum of their radii, this is:

$$\begin{array}{c}\hfill r\le R_1+R_2,\end{array}$$(2)

with r being the distance between both stars center of mass, and R₁ and R₂ are the radii of the stars. Given that the stars are physically in contact at this point, we also call this a stellar collision. Additionally, since each collision in our runs led to a merger, we then consider the terms merger and collision to be synonymous. Once this condition was fulfilled, we merged both stars replacing them by a new star whose total mass is the sum of the masses of the merging stars. Furthermore, the radius was calculated assuming that the star rapidly settles to a stable configuration where the new star has the same density as the previously merging stars, which is consistent with the calculations of Hosokawa et al. (2012) and Haemmerlé et al. (2018). Therefore the new mass and radius of the merger product were calculated as:

$$M_{\mathrm{new}}=M_1+M_2,$$(3)

$$R_{\mathrm{new}}=R_1{\left(\frac{M_1+M_2}{M_1}\right)}^{1/3}.$$(4)

3. Results

First we present the results of the simulations that do not include stellar mergers and describe the effects introduced by the addition of an external potential and the variation of the potential mass. Then, we present the results of our simulations that include stellar mergers and describe how their evolution and the formation of massive stars due to the runaway merger process are affected by the external potential.

3.1. Simulations without mergers

Here, we describe the main effects of including an external potential on the evolution of our star cluster models in the absence of stellar mergers. Due to the similarities between the runs, here, we only present the results for simulations number 5 and 8 listed in Table 1 with the setup described in Sect. 2.1. However, a description for the remaining simulations listed in Table 1 is included in Appendix B.

The evolution of clusters with M_ext = 0 × M_stars is well understood and consistent with the results of Spitzer (1987), where the clusters evolve toward core collapse as a result of two-body relaxation. The core then experiences core-oscillations, contracting and expanding again due to the formation of hard binaries that act as an energy source in the core. The core collapse is expected to occur between 10 and 20 half mass relaxation times t_rh that we calculated as (Spitzer 1987):

$$\begin{array}{c}\hfill t_{\mathrm{rh}}=0.138\frac{N}{ln(\gamma N)}{t}_{\mathrm{cross}},\end{array}$$(5)

with γ = 0.4 because we used equal mass stars and t_cross was calculated with Eq. (1) setting M_tot = M_stars. The half mass relaxation time for the cluster with N = 10 000 stars and M_ext = 0 × M_stars is t_rh = 166.38 t_cross and core collapse occurs at 2306 t_cross which is 14 t_rh. We found the core collapse time by visual inspection and looking for the first drop and subsequent rise in the 10% Lagrangian radius as marked by the vertical green line in Fig. 1.

Fig. 1. Evolution of a cluster with N = 10 000 stars with a total stellar mass of Mstars = 104 M⊙ and virial radius Rv = 1.0 pc. Top panel: fraction of ejections from the cluster and bottom panel: Lagrangian radii of the cluster. The vertical green line marks the moment of core collapse.

At the moment of core collapse, the 10% Lagrangian radius drops to 0.087 pc, yielding a mean density of 3.6 × 10⁵ M_⊙ pc⁻³. After this phase, the entire cluster begins to expand and the onset of ejections takes place with a total number of ejections of about 1500 until a time of 10 000 t_cross. A star is considered to be ejected from the cluster if its distance to the center of mass of the system is > 20 R_v and its kinetic energy is higher than its potential energy at this location.

When we included an external potential in our simulations, we note first that, since the clusters start virialized, the velocities of the stars are higher compared to clusters without the background potential. This compensates for the extra gravitational force. Then we had to calculate the crossing time of the cluster using Eq. (1) with M_tot = M_stars + M_ext.

We find that all of the clusters behave in a similar way. All of them evolve toward core collapse while there is little expansion of the outer parts (> 30% Lagrangian radii) compared to clusters without an external potential (see Fig. 2).

Fig. 2. Evolution of a cluster with N = 10 000 stars with a total stellar mass of Mstars = 104 M⊙ and a virial radius of Rv = 1.0 pc in the center of an external analytic potential with Mext = 104 M⊙. Top panel: fraction of ejections from the cluster and bottom panel: Lagrangian radii of the cluster. The vertical green line marks the moment of core collapse.

For the case when the mass of the external potential is comparable to the total mass (M_ext = M_stars), the core collapse is delayed significantly, which occurs now at 51 889 t_cross (see green vertical line in bottom panel of Fig. 2). This is 23 times (in units of the cluster crossing time) later than the core collapse time for the cluster without a background potential. The 10% Lagrangian radius at this moment drops to 0.073 pc, which yields a mean density of 6.0 × 10⁵ M_⊙ pc⁻³. For this simulation, the time is not long enough to see the contraction expansion phase that follows after core collapse. Finally, in this run, only one star is ejected from the cluster due to the higher escape velocity compared to the cluster without the external potential.

In order to understand the delay in the core collapse time when including a background potential, we define the parameter q = M_ext/M_stars as the ratio between the mass of the external potential and the total mass in stars, then we find that the core collapse time scales as t_cc ∝ (1 + q)⁴t_cross as presented in Fig. 3. If we assume, as found by Spitzer (1987) that the core collapse time is proportional to the relaxation time, then according to our data the relaxation time is proportional to t_rh ∝ (1 + q)⁴, which is consistent with Eq. (A.12) for which we derived the relaxation time of a cluster in the center of an analytic external potential (see Appendix A):

$$\begin{array}{c}\hfill t_{\mathrm{rh}}=0.138\frac{N{(1+q)}^4}{log(\gamma N)}{t}_{\mathrm{cross}},\end{array}$$(6)

Fig. 3. Core collapse time tcc divided by the crossing time tcross for clusters with N = 1000 (blue circles) and N = 10 000 (red circles) as a function of q = Mext/Mstars. We also plotted two relations that scale as (1 + q)4 which fit our data and therefore suggest that if the core collapse time is proportional to the relaxation time of the cluster, then the relaxation time scales as (1 + q)4, which is the same relation that we derive in Appendix A.

where the crossing time must be calculated from Eq. (1) with q = M_ext/M_stars (see the Appendix in Leigh et al. 2013a and Leigh et al. 2014 for similar derivations but slightly modified to consider different cases).

A second important effect that we note is on the binary population. Our results show that in increasing the mass of the background potential, the number of binary systems that form decreases dramatically as shown in Fig. 4 due to the increased mean kinetic energy of the stars. Consequently, we also find that these binaries tend to be more tightly bound (see Fig. 5) given that they form at the hard-soft boundary, which is smaller for higher velocity dispersion (Leigh et al. 2015). This will certainly have an impact on the growth of a massive star through stellar mergers given that collisions are more likely to occur in binary systems.

Fig. 4. Number of binary systems as a function of time and for different values for the mass of the external potential Mext for simulations with N = 1000 and up to 8000 crossing times. The vertical lines mark the moments of core collapse and coincides with the onset of the formation of binary systems.

Fig. 5. Distribution of the semi-major axis of the binary population for the set of simulations listed in Table 1. The data clearly shows that increasing the mass of the external potential the formed binaries tend to be more tightly bound. Despite the fact that a low mass of the external potential (Mext <  1.0 ×  Mstars) may promote the formation of wider binaries, these systems tend to be short-lived due to the higher velocities of all the stars. This figure is only intended to show the distribution of the semi-major axis and not the total number of binaries that formed.

3.2. Simulations that include stellar mergers

In the following, we present the results that we obtained when we included stellar mergers in the evolution of star clusters with and without an external Plummer potential. Additionally, we describe the effects introduced by this background potential on the formation of massive merger products. These simulations are listed in Table 2.

The first effect that we note is a delay in the runaway growth of the central object, which is explained by the delay in the core collapse time due to a larger relaxation time when increasing the mass of the potential, as can be seen in Eq. (6). Moreover with this external force, binary systems are harder to form, as shown in Fig. 4 and these binaries tend to be more compact.

Here we describe the time evolution for two clusters with N = 1000 stars, with and without an external potential, which correspond to simulations number 17 and 45 in Table 2 and we also present their evolution in Figs. 6 and 7, respectively. These simulations serve as examples of the typical evolution of the runs that include stellar mergers, given that all the simulated clusters follow very similar patterns.

Fig. 6. Evolution of a cluster with N = 1000 stars, total mass Mstars = 104 M⊙, Rv = 0.14 pc, and Rstar = 100 R⊙. Uppermost panel: number of mergers as a function of time, the fraction of binary, hyperbolic, and chain mergers in the second panel along with the fraction of stars that escape from the cluster. Third panel: 10, 50 and 90% Lagrangian radii, and bottom panel: mass growth of the most massive object that formed through mergers in the cluster.

Fig. 7. Evolution of a cluster with N = 1000 stars, stellar mass Mstars = 104 M⊙, Rv = 0.14 pc, and Rstar = 100 R⊙ in the center of an external potential with Mext = 104 M⊙. Uppermost panel: number of mergers as a function of time, the fraction of binary, hyperbolic, and chain mergers in the second panel, along with the fraction of stars that escape from the cluster. Third panel: 10, 50, and 90% Lagrangian radii, and bottom panel: mass growth of the most massive object that formed through mergers in the cluster. The drop in the 50% Lagrangian radius toward the end of the simulation suggests that the central object may grow up to 5000 M⊙.

In our runs, we identify three types of mergers: namely hyperbolic mergers, which are mergers that occur between stars that are not gravitationally bound; binary mergers, which occur between stars that are part of a binary system; and chain mergers, which occur between stars that are part of a higher order system, for example, triples. Before the onset of runaway mergers occurs in all of the clusters, there is a nearly constant merger rate for hyperbolic mergers which do not produce a single massive star, but instead several less massive stars that eventually sink to the cluster center, thus finally contributing to the growth of the most massive star. We note that the hyperbolic merger rate in clusters that do not include an external potential is only maintained for short periods of time (see Fig. 6); whereas in clusters that include the extra force, the hyperbolic merger rate is maintained until the end of the simulation (see Fig. 7). Chain mergers are very rare with only a handful having been identified among all of our simulations.

During the period before the rapid growth of the most massive star, relaxation processes drive an energy flow from the central parts of the cluster (inside the half-mass radius) via high velocity stars that migrate to the cluster halo and cause the expansion of the outer parts. All of the clusters evolve toward core collapse, but not all of them reach that stage because for large stellar radii R_star ≥ 100 R_⊙, the growth of the central star is fast and the drop in the 10% Lagrangian radius that we see in Fig. 6 occurs because the central object has already gathered 10% of the cluster mass, and thus we see the position of the central star rather than a collapsed core. The onset of growth for the central star is also marked by a rapid increase in the merger rate of stars in binary systems despite the presence of an external potential.

An important difference that we note when including the background potential is on the expansion of the clusters and the fraction of escaping stars. In this sense we found that the external potential prevents the evaporation of the cluster (see panel 2 of Fig. 7) and keeps it compact so the mergers can continue for longer periods of time, potentially gathering up to 50% or more of the cluster mass in the central object, as the third panel of Fig. 7 suggests.

3.3. Formation and growth of the most massive star

One important goal of this research is to address the effect of the external potential on the formation of a massive star by stellar mergers. Thus, we present the results regarding the mass growth of the most massive star in our clusters, along with a description of the effects introduced by the background potential.

The first important effect is a delay in the formation of the most massive star due to the increased stellar velocity, which in turn causes relaxation processes to be slower. This then leads to a reduction in the number of binary systems that form, at a given time, if the mass of the external potential is comparable to the total stellar mass. On the other hand, the external potential can also favor the mass growth of the central star by preventing the evaporation and expansion of the cluster.

In the following, we present a model that we used to estimate the mass of the central object at different times for both a cluster with and without an external potential. For this, we followed the same method described in Reinoso et al. (2018).

In the lower panel of Fig. 8, we first present an example of the mass growth of the most massive star, which undergoes a very rapid growth at around 1000–1500 t_cross. This coincides with the peak in the number of mergers with the central object as shown in the top panel of the same figure.

Fig. 8. Example of a Gaussian fit to the number of mergers with the central star for simulation number 45 listed in Table 2. Top panel: number of mergers with the most massive star as a function of time (black solid line) along with the Gaussian function that we used to estimate the total number of mergers at different times. The parameters tdelay and tduration are also shown in the figure. Bottom panel: mass growth of the central star whose rapid growth coincides with the peak in the number of mergers.

We fit the Gaussian function presented in Eq. (7) to the combined data from the six random realizations for each simulation setup presented in Table 2 in order to get an estimate for the number of mergers with the central object during the rapid growth. Thus, we define the delay time t_delay as the time at which the peak in the Gaussian occurs, and the height A of the Gaussian gives us an estimate for the number of mergers at t_delay. Additionally, we define the duration time $t_{\mathrm{duration}}=\mathit{FWHM}/\left(2\sqrt{2ln2}\right)$, with FWHM being the full width at half maximum of the Gaussian. An illustration of the fitting formula and two parameters are shown in the top panel of Fig. 8. By doing this, we were able to obtain an estimate for the moment at which there is rapid growth of the most massive star (t_delay) and an estimate for the duration of this period (t_duration), which we can later combine with A and a normalization factor D to finally obtain the number of mergers experienced by the central star using Eq. (8) and therefore an estimate for its final mass at different times. Furthermore we determined D by comparing our model to the results from our simulations. We also show how well our model reproduces the data of our simulations in Figs. 9 and 10. The total number of mergers at time t can thus be calculated as:

$$\begin{array}{c}\hfill N_{\mathrm{merger}}(t)=Aexp(-\frac{{(t-t_{\mathrm{delay}})}^2}{2t_{\mathrm{duration}}^2}),\end{array}$$(7)

Fig. 9. Mean mass of the most massive star normalized by its initial mass, as a function of N and Rstar for clusters without an external potential evolved until 2000 tcross = 15.6 Myr. The symbols represent the mean values obtained from six simulations and the solid lines show the results of the calculations with our model presented in Eq. (8) using t = 15.6 Myr and the parameters obtained from Eqs. (9)–(12).

Fig. 10. Mean mass of the most massive star normalized by its initial mass, as a function of N and Rstar for clusters without an external potential evolved until 3000 tcross = 16.5 Myr, except for the point at N = 5000 and Rstar = 20 R⊙, which are simulations that were evolved until 4000 tcross = 22 Myr. The symbols represent the mean values obtained from six simulations and the solid lines show the results of the calculations with our model presented in Eq. (8) using t = 16.5 Myr (t = 22 Myr for N = 5000 and Rstar <  50 R⊙) and the parameters obtained from Eqs. (13)–(16).

and the total number of mergers until a time t_end can be calculated as:

$$\begin{array}{c}\hfill {\displaystyle N_{\mathrm{merger},\mathrm{total}}=AD{\int }_0^{t_{\mathrm{end}}}exp(-\frac{{(t-t_{\mathrm{delay}})}^2)}{2t_{\mathrm{duration}}^2})\mathrm{d}t.}\end{array}$$(8)

Although the data is not a perfect Gaussian (see Fig. 8), we only need an estimate for the delay time t_delay, the duration time t_duration, and an estimate for the number of collisions A at t_delay. Following this procedure, we first combined the data of the six random realizations per each of the configurations listed in Table 2 and then applied the Gaussian fit of Eq. (7) to find the values of A, t_delay, and t_duration that we present in Figs. 11 and 12 for clusters with M_ext/M_stars = 0.0 (simulations 1–28 in Table 2) and in Fig. 13 for clusters with M_ext/M_stars = 1.0 (simulations 29–56 in Table 2). Then we fit these data using an implementation of the nonlinear least-squares Marquardt-Levenberg algorithm in gnuplot and obtained Eqs. (9)–(11) for clusters without an external potential and Eqs. (13)–(15) for clusters in a background potential. Using these equations, we were able to compute N_{merger, total} assuming first D = 1 and we compared the calculated values to the real values from the simulations to adjust D and reproduce the results correctly. By doing this, we obtain Eqs. (12) and (16):

$$\begin{array}{cc}& {log}_{10}(A)=[0.104{log}_{10}(R_{\mathrm{star}})+0.562]{log}_{10}(N)\hfill \\ \hfill & +0.358{log}_{10}(R_{\mathrm{star}})-1.011,\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill & {log}_{10}(t_{\mathrm{delay}}[t_{\mathrm{cross}}])=[-0.246{log}_{10}(N)+0.443]{log}_{10}(R_{\mathrm{star}})\hfill \\ \hfill & +0.954{log}_{10}(N)+0.077,\hfill \end{array}$$(10)

$${log}_{10}(t_{\mathrm{duration}}[t_{\mathrm{cross}}])=0.430{log}_{10}(N)+0.580,$$(11)

$${log}_{10}(D)=-0.227{log}_{10}(R_{\mathrm{star}})-1.791,$$(12)

Fig. 11. Parameters of the Gaussian fit to the number of mergers with the central star for simulations of clusters without an external potential. The data points were obtained after the fit to the combined data of six simulations. Upper panel: values of the normalization factor A for different values of N and Rstar along with our model from Eq. (9). These data points show that the number of mergers increases with both N and Rstar as expected. Bottom panel: different values of tdelay in addition to the different values of N and Rstar along with the fitting function presented in Eq. (10). These points show that tdelay decreases with Rstar, but it increases with N.

Fig. 12. tduration as a function of the number of stars N for clusters without an external potential. We found no clear correlations between tduration and Rstar. The solid blue line shows the fitting function from Eq. (11).

Fig. 13. Parameters of the Gaussian fit to the number of mergers with the central star for simulations of clusters with an external potential. The data points were obtained after the fit to the combined data of six simulations. Upper panel: values of the normalization factor Aext for different values of N and Rstar along with our model from Eq. (8). These data points show that the number of mergers increases with both N and Rstar as expected. Middle panel: different values of tdelay, ext in addition to the different values of N and Rstar along with the fitting function presented in Eq. (14). These points show that tdelay, ext decreases with Rstar, but it increases with N. Additionally, these values are larger than the values of tdelay for clusters without an external potential (see Fig. 11). Lower panel: values of tduration, ext along with the fitting function from Eq. (15). We clearly see a positive correlation with N and a negative correlation with Rstar.

and the parameters for clusters in an external potential:

$$\begin{array}{cc}& {log}_{10}(A_{\mathrm{ext}})=[3.975\times {10}^{-5}N+0.606]{log}_{10}(R_{\mathrm{star}})\hfill \\ \hfill & +0.554{log}_{10}(N)-1.510,\hfill \end{array}$$(13)

$$\begin{array}{cc}\hfill & {log}_{10}(t_{\mathrm{delay},\mathrm{ext}}[t_{\mathrm{cross},\mathrm{ext}}])=[-0.006{log}_{10}(N)+0.240]{log}_{10}(R_{\mathrm{star}})\hfill \\ \hfill & +0.064{log}_{10}(N)-1.170+t_{\mathrm{delay}},\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill & {log}_{10}(t_{\mathrm{duration},\mathrm{ext}}[t_{\mathrm{cross},\mathrm{ext}}])=-0.178N^{0.119}{log}_{10}(R_{\mathrm{star}})\hfill \\ \hfill & +2.163N^{0.073}+0.151,\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill & {log}_{10}(D_{\mathrm{ext}})=0.053{log}_{10}(R_{\mathrm{star}})\hfill \\ \hfill & -1.879\times {10}^{6.9N/{10}^6}-0.778.\hfill \end{array}$$(16)

It is important to note that in Eqs. (10) and (11), the delay time t_delay and duration time t_duration are expressed in units of the crossing time of the cluster defined in Eq. (1) and consequently when calculating the total number of mergers until a time t_end using Eq. (8), this time must also be expressed in units of the crossing time t_cross. The same principle applies for Eqs. (14) and (15). That is, when calculating the number of mergers until a time t_end using Eq. (8), this time must be expressed in units of the cluster crossing time that have been modified by the mass of the external potential t_cross, ext, which is defined in Eq. (1).

The delay time for clusters in a background potential is defined in Eq. (14) as t_delay, ext = t_delay + t_{delay, Mext = Mstars}, which is the delay time for clusters without an external potential t_delay plus an additional term t_{delay, Mext = Mstars} = [ − 0.006log₁₀(N)+0.240]log₁₀(R_star)+0.064log₁₀(N) − 1.170 that, in principle, may depend on the mass of the potential. However we do not examine this potential dependence here.

To help in the reading of the equations, we define a general equation of the form:

$$\begin{array}{c}\hfill {log}_{10}(x)=\alpha {log}_{10}(R_{\mathrm{star}})+\beta {log}_{10}(N)+\gamma \end{array}$$(17)

and present the parameters α, β, and γ in Table 3. They can be used to recover Eqs. (9)–(16).

Now we can estimate the number of mergers with the central star and thus the mass of this star up to a time t for clusters with different initial conditions, including the presence of an external potential that is comparable to the total stellar mass in the cluster. This can be used as a very simplified model of a star cluster, which still contains gas during or shortly after the process of star formation or a nuclear star cluster. In Fig. 14, we show the expected fraction f_merger of the stars that merge with a single central object in a cluster without an external potential. In Fig. 15, this is shown for clusters in an external potential, for a broad combination of N and R_star after 1, 2, and 10 Myr assuming a cluster with M_stars = 10⁴ M_⊙, a virial radius of R_v = 0.14 pc, and a mass of the external potential M_ext = M_stars for the second case. We used Eq. (8) for these calculations along with the parameters from Eqs. (9)–(12) to obtain the results presented in Fig. 14, and Eqs. (13)–(16) to obtain the results presented in Fig. 15. Our model gives the total number of stars N_{merger, total} that merge with the central object in a given time interval, and thus assuming that all stars are equal during the time at which most of the collisions occur, we can estimate the mass of the central object as $M_{\max }=N_{\mathrm{merger},\mathrm{total}}\overline{m}=N_{\mathrm{merger},\mathrm{total}}{M}_{\mathrm{stars}}/N=f_{\mathrm{merger}}{M}_{\mathrm{stars}}$. Our results indicate that when collisions are maintained for short periods of time only, that is, 1 Myr, then the clusters without an external potential form more massive central stars than clusters with an external potential (see left panels of Figs. 14 and 15). This is due to the fact that the external potential dramatically increases the relaxation time of the cluster, as shown in Eq. (6). In fact the relaxation time depends on the ratio of external potential mass to stellar mass q = M_ext/M_stars as t_relax ∝ q⁴.

Fig. 14. Fraction of stars fmerger that merge with the central star in a cluster with a total stellar mass of Mstars = 104 M⊙ and a virial radius of Rv = 0.14 pc after 1 Myr (left panel), 2 Myr (middle panel), and 10 Myr (right panel) for different combinations of the number of stars N and the stellar radius Rstar calculated using our model presented in Eq. (8). We can see that when the collisions are maintained for short periods of time, that is, 1 Myr, a massive object with ∼fmergerMstars = 0.15 × 104 M⊙ = 1500 M⊙ is formed in the cluster center only when N is small (given that tdelay becomes larger than the integration time for large N) or large Rstar. If the time limit is long, that is, 10 Myr, most of the clusters will form stars with ∼0.3 Mstars = 3000 M⊙.

Fig. 15. Fraction of stars fmerger that merge with the central star in a cluster with a total stellar mass of Mstars = 104 M⊙, a virial radius of Rv = 0.14 pc, and an external potential with a mass of Mext = Mstars after 1 Myr (left panel), 2 Myr (middle panel), and 10 Myr (right panel) for different combinations of the number of stars N and the stellar radius Rstar calculated using our model presented in Eq. (8). We can see that when the collisions are maintained for short periods of time, that is, 1 Myr, an object that is not so massive with ∼fmergerMstars = 0.1 × 104 M⊙ = 1000 M⊙ is formed in the cluster center only when N is very small (given that tdelay becomes larger than the integration time for large N) or when Rstar is large. If the time limit is 2 Myr (middle panel), we expect that clusters with a small N experience core collapse and form a massive central star. However, if the number of stars is large and Rstar is also large, given the higher velocity of the stars in those clusters compared to clusters without the external potential, there is a larger probability for hyperbolic mergers to occur before core collapse, which we see on the right-hand side of middle panel. If the time limit is long, that is, 10 Myr, most of the clusters form stars with ∼0.3−0.4 Mstars = 3000−4000 M⊙ and even 0.7 Mstars = 7000 M⊙ if N is large enough.

As the time limit increases, we see that in both clusters with and without an external potential a more massive star emerges in the center. However, even more massive stars form in the cluster without an external potential for most values of N. However, there is an important difference in clusters that include a background potential, which is that the potential keeps the cluster compact as previously found by Leigh et al. (2014), and thus most of the cluster mass is able to eventually sink to the center at later times. This effect is visible in the left part of the middle panel in Figs. 14 and 15. We also note that for large N and R_star for clusters in a background potential and before core collapse, we may also expect mergers with the central star, not binary mergers but hyperbolic, given that in those cases the cross section for collisions is very large and the probability of hyperbolic collisions increases with stellar velocity, that is, if there is an external potential. Finally when the time limit is very long, that is, 10 Myr, we may expect final masses for the central object in the order of 0.3–0.4 M_stars for most clusters without an external potential and 0.4–0.7 for clusters in a background potential because the compactness of the cluster is maintained and collisions still occur at an approximately constant rate after the stage of runaway growth of the central star.

4. Implications for primordial clusters

In the following, we explore the implications of our results with respect to primordial star clusters, including both embedded and gas free. Here we particularly distinguish the case of standard Pop. III clusters as expected in a typical minihalo with about 10⁶ M_⊙.

4.1. Standard Pop. III clusters (minihalos)

For a typical Pop. III star cluster, we assume a mass of 1000 M_⊙, which is consistent with a baryon fraction of about 10% and a star formation efficiency on the order of 1% in a 10⁶ M_⊙ minihalo. We adopt a virial radius for the cluster on the order of 0.1 pc, which is consistent with the results from simulations and semi-analytic models (Clark et al. 2011a,b; Greif et al. 2011, 2012; Latif et al. 2013b; Latif & Schleicher 2015). We use a stellar radius of 100 R_⊙, which is characteristic of primordial protostars with accretion rates on the order of 10⁻³ M_⊙ yr⁻¹ (Hosokawa et al. 2012). The crossing time of the cluster then is 0.015 Myr. The number of stars that can be expected in such a cluster is uncertain, but here we adopt an estimate of about N = 100. Thus, a stellar mass of m_star = 10 M_⊙. Using the relations that we found, we expect two collisions to occur within 1 Myr which correspond to a final mass of 30 M_⊙. If we assume now a lifetime of 10 Myr, we can expect a total of six collisions with the central star which correspond to a final mass of 70 M_⊙.

We consider now an embedded and accreting Pop. III star cluster with a stellar mass of 1000 M_⊙, a gas mass of 1000 M_⊙, a virial radius of 0.1 pc, and a mean stellar radius of 100 R_⊙. The crossing time of the cluster, according to Eq. (1), corresponds to 0.007 Myr. We also adopt a number of stars N = 100. Then our model predicts that no collisions occur within 1 Myr, but for a lifetime of 10 Myr we expect a total of 15 collisions and hence a central star with 160 M_⊙. We expect the lifetime of a massive primordial star to be in between this range, depending on the precise mass, the amount of rotation, and the effects that the collisions may have on the stellar evolution (Maeder & Meynet 2012). We therefore find that a moderate enhancement can be achieved within a normal cluster. We note that the values given here are the expected mean number of mergers. Individual clusters can deviate from these, both towards lower and higher fractions of mergers, potentially including clusters with zero mergers. This is especially the case for typical Pop. III clusters where the number of collisions can be expected to be comparable to the mean value because it is a chaotic collisional dynamical process.

4.2. Massive primordial clusters (atomic cooling halos)

As a next step, we address now the potential impact of collisions in a more massive atomic cooling halo with a total mass of 10⁸ M_⊙. Under the right conditions and in particular if the cooling on larger scales is regulated by atomic hydrogen (Latif et al. 2014), a rather massive cluster of 10⁴ M_⊙ can form, which is then exposed to larger accretion rates on the order of 10⁻¹ M_⊙ yr⁻¹. We assume that the cluster consists of an initial number N = 1000 of stars, the virial radius is R_v = 0.14 pc, and that the stellar radii are somewhat enhanced compared to the standard Pop. III cluster due to the higher accretion rates, with a typical radius of about 300 R_⊙. The crossing time of the cluster according to Eq. (1) with M_tot = 10⁴ M_⊙ is then 0.0078 Myr. Using our derived model we expect about 51 collisions in 1 Myr, and about 168 within 10 Myr with a single central star. We again expect the realistic lifetime of the resulting massive star to be in between these extreme cases. In the case of an atomic cooling halo, we thus conclude that a considerable enhancement is possible as a result of stellar mergers. If we take the mean stellar mass, then the expected masses are 520–1690 M_⊙ after 1 and 10 Myr, respectively.

If we now consider the same cluster configuration but during the embedded phase, assuming a total mass in gas of 10⁴ M_⊙, the crossing time of the cluster according to Eq. (1) with q = M_ext/M_stars = 2 is then 0.0039 Myr, and we expect a total of 16 collisions to occur within 1 Myr and a total of 199 mergers within 10 Myr. Again, by taking the mean stellar, mass this corresponds to a mass of 170–2000 M_⊙ within 1 and 10 Myr, respectively. Also the reported values here correspond to a mean, and there can be deviations to lower and higher merger fractions. Regardless, we expect the number of collisions and hence the final mass to remain within the same order of magnitude.

5. Discussion

We modeled the runaway growth of stars through mergers in the center of dense star clusters by including an analytic external potential in our N-body simulations. Our model relies on computing four parameters that depend on the number of stars N and their stellar radii R_star. These parameters are the delay time t_delay, the duration time for the collision process t_duration, and two normalization factors A and D. We then used these four parameters to integrate a Gaussian function and estimate the mass of the central object up to a time t. Although a Gaussian fit is a good tool to estimate t_delay, t_duration, and A (see Fig. 8), there is a deviation from this function when we include an external potential. In fact, when we compare the mass growth of the central object in clusters with and without an external potential (see, e.g. Figs. 6 and 7), we see a delay for clusters in an external potential and an initially less dramatic growth until the onset of mergers of stars in binary systems. This growth continues after the runaway growth and is seen as tails of a Gaussian in a plot of N_mergers versus time (see top panel of Figs. 9 and 10). This later growth is associated to hyperbolic mergers, that is, mergers of stars that are not bound by gravity, and this type of merger becomes important when including an external potential because the cluster remains compact for longer periods of times (Leigh et al. 2014). The onset of collisions is delayed in simulations with an external potential relative to simulations without it.

We also have found a formula for the delay time when we include an external potential t_delay, ext = t_delay + t_{delay, Mext = Mstars}, which is basically the same expression for clusters without the external potential plus an additional term t_{delay, Mext = Mstars} presented in Eq. (14). We present this parameter in this way so that we can investigate in the future if this can be used in a more general way as a function of the mass of the external field M_ext and if it is related to a longer relaxation time by means of its dependence on M_ext as suggested by Eq. (6).

The model derived here can be used as a first approximation to obtain the number of stellar mergers in embedded star clusters or to understand the effects of an external potential on the formation of massive merger products, which might also be important for the modelling of nuclear star clusters. Future research in this field will employ hydrodynamic modeling of the gas, which may prevent a large delay in the onset of runaway growth for the central star due to gas accretion and dissipative effects via star-gas or binary-gas interactions. These simulations also need to include realistic mass-radius relations for accreting protostars when aiming for a realistic modeling of the first star clusters in the Universe. By performing these suites of N-body simulations we have covered a large parameter space, and this allows us to provide hints as to which part of this parameter space the future, more sophisticated simulations should be focused on and the subsequent physics that becomes relevant and often even dominant in these regimes.

Our simulations including mergers adopt a virial radius of R_v = 0.14 pc and a total stellar mass of M_stars = 10⁴ M_⊙, but our results can be rescaled for different sizes and masses by means of the crossing time. In this work, we only consider for the initial conditions equal mass and equal radii stars but allow for the mass and radius to vary due to the mergers according to Eqs. (3) and (4), and in this sense our simulations represent an effective model where an average stellar radius is adopted over the period of time considered. This effective stellar radius should correspond to the typical radius when the majority of mergers are expected to occur.

Our model still lacks a more realistic stellar population, namely an initial mass function (IMF), which would naturally lead to mass segregation and subsequent mergers of the most massive stars in the cluster center. Moreover, we have not included the dissipative effects expected for star-gas interactions in our simulations that include an external potential, which could lead to the formation of more binary systems and hence more mergers (Leigh et al. 2014). In this sense, we do not expect a significant variation in the final mass of the most massive stars if the total mass in the most massive stars is comparable to the final mass obtained in our models, around 30–2000 M_⊙. Additionally, including an IMF, three-body encounters and two-body relaxation would cause the ejection of the smallest stars in the cluster, which could be an important source of stellar relics from the first star clusters in the early Universe. Another important way to proceed is a better modeling of the stellar mergers, in particular the mass loss during the stellar collisions which may be up to 25% of the total mass (Gaburov et al. 2010) and the collision product may even receive a kick velocity > 10 km s⁻¹. More recent work on the effects of mass loss and its effect on runaway growth of a central star in a cluster shows that considering 5% mass loss in every collision produces objects which are 20–40% less massive at the end of the runaway growth compared to models without mass loss (Glebbeek et al. 2009; Alister Seguel et al. 2020). Thus regarding this, the results presented here should be taken as an upper limit on the mass of the merger product.

6. Conclusions

In this study, we have performed a set of 374 N-body simulations of dense, virialized star clusters and clusters in the center of an external potential including stellar mergers in order to understand the effects of the background force on the formation of massive stars and derive a model to estimate the mass enhancement of these objects in embedded star clusters. We find that the presence of an external potential delays the overall evolution of the star cluster and the formation of a massive central star through runaway collisions due to the increased kinetic energy of the stars, which in turn increases the relaxation time. However, the merger products become more massive (if the collision process is maintained for a long time) given that these clusters expand less than clusters without an external potential and so more stars are able to merge with the central star even after the process of runaway growth.

We also find that the increased velocity dispersion for star clusters in an external potential boosts the number of hyperbolic mergers , that is, mergers of stars that are not part of a binary or triple system, with up to 50–60% of the total number of mergers being hyperbolic mergers. Whereas in clusters without an external potential, this percentage is around 30%. This is due to the fact that in presence of a background potential, the clusters remain more compact (Leigh et al. 2013a, 2014) and hyperbolic mergers still may occur. Stellar ejections are highly suppressed in clusters in the center of an external potential.

We find a set of equations that can be used to estimate the mass of the merger product at different times for both clusters with and without an external potential, and we present an example of such calculations in Figs. 14, 15, and in Sect. 4 proving that if the process is not interrupted for long periods of time (≥10 Myr), the external potential enhances the mass of the central object by a factor of ∼2. However, if the process is interrupted at early times (∼1 Myr), the clusters in a background potential produce objects only half as massive compared to objects formed in clusters without an external potential. When applied to Pop. III star clusters we find, for standard clusters formed in a minihalo, a moderate enhancement for the mass of the most massive star in the range of 10–160 M_⊙ within 1 and 10 Myr in embedded clusters. Whereas for a massive Pop. III star cluster that formed in an atomic cooling halo, we find that the mass of the most massive star lies in the range of 170–2000 M_⊙ in embedded clusters within 1 and 10 Myr.

Acknowledgments

We thank the anonymous referee for the useful comments. BR thanks Sverre Aarseth for his help with the code NBODY6. BR also thanks Conicyt for funding through (CONICYT-PFCHA/MagísterNacional/2017-22171385), (CONICYT-PFCHA/Doctorado acuerdo bilateral DAAD/62180013), financial support from DAAD (funding program number 57451854). DRGS acknowledges support of Fondecyt Regular 1161247 and through Conicyt PIA ACT172033. BR and DRGS acknowledge financial support through Anillo (project number ACT172033). MF acknowledges funding through Fondecyt Regular No.1180291. NWCL gratefully acknowledges the support of a Fondecyt Iniciacion #11180005.

References

Appendix A: Modification of the half-mass relaxation time-scale

We have seen that when we include an external potential in a star cluster, the overall evolution seems to be delayed, in particular the core collapse times-scale rapidly increases with the mass of the external potential. Until now, we have used the usual half-mass relaxation time-scale defined in Eq. (5) (Spitzer 1987). However, this is probably not adequate when we include an external potential because in that case, considering a virialized cluster, the root mean square (rms) velocity of the stars is higher due to the extra force exerted by the potential and this certainly modifies the time-scale on which stellar encounters are going to modify the higher velocity of the stars. In order to account for this effect, here, we derive a modified relaxation time-scale with the aid of a new parameter q = M_ext/M_stars. We begin with the usual derivation of the relaxation time assuming virial equilibrium, a similar derivation can be found in Binney & Tremaine (1987):

$$\begin{array}{c}\hfill 2T+U=0.\end{array}$$(A.1)

We note that $T=\frac{1}{2}{M}_{\mathrm{stars}}⟨v^2⟩$ is the kinetic energy and $|U|=\frac{\mathit{GM}}{2R_{\mathrm{v}}}$ is the potential energy. In this case, M_stars is the total mass in stars, ⟨v²⟩ is the mean velocity squared of the stars, M = M_stars + M_ext is the total mass, that is, the mass in stars plus the mass of the external potential M_ext, and R_v is the virial radius of the cluster, which in this case is the same virial radius for the external potential. Now we write the total mass as a function of q:

$$\begin{array}{cc}& M=M_{\mathrm{stars}}+M_{\mathrm{ext}},\hfill \\ \hfill & M=M_{\mathrm{stars}}(1+\frac{M_{\mathrm{ext}}}{M_{\mathrm{stars}}}),\hfill \\ \hfill & M=M_{\mathrm{stars}}(1+q).\hfill \end{array}$$(A.2)

We begin now with a derivation for the change of the velocity of a star due to an encounter with another star. The force that one star feels due to the other star is (see Fig. A.1):

$$\begin{array}{c}\hfill F=\frac{\mathit{GMm}}{d^2}=\frac{\mathit{GMm}}{b^2+V^2{t}^2},\end{array}$$(A.3)

Fig. A.1. Trajectory of a star with mass M that passes close to another star of mass m, which causes a small deflection on the velocity vector V.

where M is the mass of one star and m is the mass of the other star, and d is the impact parameter or distance of closest approach. Now we assume that the velocity change ΔV of the star with mass M in the direction parallel to V is small, and only the perpendicular component of the velocity is changed due to the perpendicular force:

$$\begin{array}{cc}& F_{\perp }=Fsin\varphi =F\times \frac{b}{d}=\frac{\mathit{GMmb}}{{(b^2+V^2{t}^2)}^{3/2}},\hfill \end{array}$$(A.4)

$$\begin{array}{cc}& F_{\perp }=M\frac{\mathrm{d}{V}_{\perp }}{\mathrm{d}t}.\hfill \end{array}$$(A.5)

Next we want to know the time at which the perpendicular velocity changes by an amount ΔV ∼ V, which is the relaxation time-scale. Therefore, we integrate Eq. (A.5) so the final perpendicular velocity is:

$$\begin{array}{cc}& {\displaystyle \mathrm{\Delta }V={\int }_{-\infty }^{\infty }\frac{\mathrm{d}{V}_{\perp }}{\mathrm{d}t}\mathrm{d}t,}\hfill \\ \hfill & {\displaystyle ={\int }_{-\infty }^{\infty }\frac{F_{\perp }}{M}\mathrm{d}t,}\hfill \\ \hfill & {\displaystyle ={\int }_{-\infty }^{\infty }\frac{\mathit{Gmb}}{{(b^2+V^2{t}^2)}^{3/2}}\mathrm{d}t,}\hfill \\ \hfill & =\frac{2Gm}{\mathit{bV}}.\hfill \end{array}$$

As the star receives many deflections in different directions, we are interested in the mean value of the squares of the velocity kicks $\langle \mathrm{\Delta }{V}_{\perp }^2\rangle$ that can be found integrating all the small deflections as:

$$\begin{array}{c}\hfill {\displaystyle ⟨\mathrm{\Delta }{V}_{\perp }^2⟩={\int }_{b_{\min }}^{b_{\max }}{\left(\frac{2Gm}{\mathit{bV}}\right)}^{2}dN,}\end{array}$$

where dN is the expected number of encounters that occurs in a time t between impact parameters b and b + db for a star with typical velocity V, this is:

$$\begin{array}{c}\hfill dN=n\times Vt\times 2\pi bdb,\end{array}$$

where n is the number density of stars, then,

$$\begin{array}{cc}& {\displaystyle ⟨\mathrm{\Delta }{V}^2⟩={\int }_{b_{\min }}^{b_{\max }}nVt{\left(\frac{2Gm}{\mathit{bV}}\right)}^{2}2\pi bdb,}\hfill \\ \hfill & {\displaystyle =\frac{8\pi G^2{m}^{2}nt}{V}{\int }_{b_{\min }}^{b_{\max }}\frac{\mathit{db}}{b},}\hfill \\ \hfill & =\frac{8\pi G^2{m}^{2}nt}{V}ln\left(\frac{b_{\max }}{b_{\min }}\right).\hfill \end{array}$$(A.6)

After enough time, the perpendicular velocity of one star grows to its original speed, this time is the relaxation time-scale that we can derive using Eq. (A.6):

$$\begin{array}{c}\hfill t_{\mathrm{relax}}=\frac{V^3}{8\pi G^2{m}^{2}nln\left(\frac{b_{\max }}{b_{\min }}\right)}.\end{array}$$(A.7)

The term $\frac{b_{\max }}{b_{\min }}$ is often written as Λ and is the ratio of the size of the system to the “strong encounter distance” b_min = 2Gm/V², which is the distance at which an encounter with another star would result in a 90 degree deflection. This ratio $\frac{b_{\max }}{b_{\min }}$, for a system of identical N stars is found to be $\frac{b_{\max }}{b_{\min }}=\gamma N$ (Spitzer 1987), with γ = 0.4. This definition of the relaxation time-scale is useful to estimate the core collapse time, which typically occurs between 15–20 half-mass relaxation times. The problem comes when we include an analytic external potential and evolve the cluster at the center of this potential. In this case, we have an extra force acting on the stars, and therefore a virialized cluster under the influence of this external potential has a larger potential energy and also a larger kinetic energy compared to a cluster in which there is not an external potential. This should increase the relaxation time-sale as the velocity of the stars increases and the number and masses of stars is kept constant. We now derive the relaxation time for a cluster under the influence of an external potential, assuming that this potential follows the same mass distribution of the stars and with the same virial radius. In Recalling the condition for virial equilibrium from Eq. (A.1) and the total mass for the cluster from Eq. (A.2) we find that:

$$\begin{array}{cc}& 2T+U=0,\hfill \\ \hfill & 2\left(\frac{1}{2}Nm{V}^2\right)-\frac{G{M}_{\mathrm{tot}}^2}{R_{\mathrm{v}}}=0,\hfill \\ \hfill & M_{\mathrm{stars}}{V}^2=\frac{G{M}_{\mathrm{stars}}^2{(1+q)}^2}{R_{\mathrm{v}}},\hfill \\ \hfill & V=\sqrt{\frac{G{M}_{\mathrm{stars}}}{R_{\mathrm{v}}}}(1+q),\hfill \end{array}$$

where V is the rms velocity of the stars, N is the number of stars, m is the mass of a single star, M = M_stars + M_ext is the total mass, that is, the mass of the stars plus the mass of the external potential, and $q=(1+\frac{M_{\mathrm{ext}}}{M_{\mathrm{stars}}})$ and R_v is the virial radius of the system. When q = 0, we get the typical velocity for stars in a cluster which is in virial equilibrium:

$$\begin{array}{c}\hfill V=\sqrt{\frac{G{M}_{\mathrm{stars}}}{R_{\mathrm{v}}}},\end{array}$$(A.8)

but for a cluster in virial equilibrium and with a background potential, the typical velocity is modified as:

$$\begin{array}{cc}& V_{\mathrm{ext}}=\sqrt{\frac{G{M}_{\mathrm{stars}}}{R_{\mathrm{v}}}}(1+q),\hfill \\ \hfill & V_{\mathrm{ext}}=V(1+q),\hfill \end{array}$$(A.9)

where V_ext is the rms velocity of the stars in the presence of an external potential and V is the rms velocity of stars without a background potential.

The relaxation time is often expressed as a function of the crossing time t_cross of the cluster, which is simply the time it takes for a star with the typical velocity V to cross the system:

$$\begin{array}{c}\hfill t_{\mathrm{cross}}=\frac{R_{\mathrm{v}}}{V},\end{array}$$(A.10)

then, combining Eqs. (A.7) and (A.10) the relaxation time for clusters without an external potential is:

$$\begin{array}{c}\hfill t_{\mathrm{relax}}=\frac{V^4}{8\pi G^2{m}^{2}n{R}_{\mathrm{v}}ln\left(\frac{b_{\max }}{b_{\min }}\right)}{t}_{\mathrm{cross}}.\end{array}$$(A.11)

Now, if we include an external potential, q ≠ 0, and in substituting Eqs. (A.9) for (A.10), the crossing time for a cluster in an external potential is:

$$\begin{array}{cc}& t_{\mathrm{cross},\mathrm{ext}}=\frac{R_{\mathrm{v}}}{V(1+q)},\hfill \\ \hfill & t_{\mathrm{cross},\mathrm{ext}}=\frac{t_{\mathrm{cross}}}{(1+q)}.\hfill \end{array}$$(A.12)

Subsequently the relaxation time for clusters in an external potential must be:

$$\begin{array}{c}\hfill t_{\mathrm{relax},\mathrm{ext}}=\frac{V_{\mathrm{ext}}^4}{8\pi G^2{m}^{2}n{R}_{\mathrm{v}}ln\left(\frac{b_{\max }}{b_{\min }}\right)}{t}_{\mathrm{cross},\mathrm{ext}}.\\ \hfill \end{array}$$

Now we replace Eq. (A.9) for Eq. (A.13) and we find:

$$\begin{array}{c}\hfill t_{\mathrm{relax},\mathrm{ext}}=\frac{V^4{(1+q)}^4}{8\pi G^2{m}^{2}n{R}_{\mathrm{v}}ln\left(\frac{b_{\max }}{b_{\min }}\right)}{t}_{\mathrm{cross},\mathrm{ext}},\end{array}$$(A.13)

Now in recalling that from the usual definition of the relaxation time we can define the half-mass relaxation time for a cluster of N equal mass stars as (Spitzer 1987):

$$\begin{array}{c}\hfill t_{\mathrm{rh}}=0.138\frac{N}{ln(\gamma N)}{t}_{\mathrm{cross}}\end{array}$$

then, for comparison purposes, we define the half-mass relaxation time-scale for clusters in an external potential as:

$$\begin{array}{c}\hfill t_{\mathrm{rh},\mathrm{ext}}=0.138\frac{N{(1+q)}^4}{ln(\gamma N)}{t}_{\mathrm{cross},\mathrm{ext}},\end{array}$$(A.14)

with $q=\frac{M_{\mathrm{ext}}}{M_{\mathrm{stars}}}$ and γ = 0.4 for equal mass stars.

Appendix B: Evolution of clusters in a background potential

In this section, we describe in more detail the results of our simulations that do not include stellar collisions and which are listed in Table 1. We show the evolution of the clusters with and without an external potential, and we find a delay in the overall evolution when increasing the mass of the external potential.

All of the clusters evolve toward core collapse, which we found upon visual inspection by looking for the first drop and subsequent rise in the 10% Lagrangian radius. The core collapse occurs for the clusters without a background potential at 456 t_cross = 19.78 t_rh for the cluster with N = 1000 (see Fig. B.1), and for the cluster with N = 10 000 at 2306 t_cross = 13.85 t_rh (see Fig. B.1).

Fig. B.1. Evolution of a cluster with N = 1000 stars, total mass Mstars = 104 M⊙ and Rv = 1.0 pc. Top panel: fraction of stars that were ejected from the cluster. Bottom panel: Lagrangian radius and the vertical green line marks the moment of core collapse. The time is presented in units of the crossing time of the cluster.

As the clusters evolve toward core collapse, the outer parts continually expand which leads, along with the onset of stellar ejections, to the evaporation of the clusters. We also note that the onset of ejections occur just after core collapse when binary systems form in the core and three body interactions lead to the ejection of stars. An ejection occurs if the distance from the center of mass of the cluster to a star is larger than 20 times the virial radius R_v and if the total energy of the star is greater than zero, that is, E_star = E_kin +|E_pot| > 0 with E_kin and E_pot being the kinetic and potential energy, respectively. These results are shown in Figs. 1 and B.1.

When we include a background potential with a low mass compared to the total mass in stars, M_ext = 0.1 × M_stars in our simulations with N = 1000, the core collapse time is 460 t_cross, which is similar to the core collapse time for the cluster without the external potential in terms of the crossing time of the cluster. However, this cluster show less expansion and only one star is ejected (see Fig. B.2).

Fig. B.2. Evolution of a cluster with N = 1000 stars, total mass Mstars = 104 M⊙ and Rv = 1.0 pc in a background potential with mass Mext = 103 M⊙. Top panel: fraction of stars that were ejected from the cluster. Bottom panel: Lagrangian radius and the vertical green line marks the moment of core collapse. The time is presented in units of the crossing time of the cluster.

Fig. B.3. Evolution of a cluster with N = 1000 stars, total mass Mstars = 104 M⊙, Rv = 1.0 pc in a background potential with mass Mext = 5 × 103 M⊙. Top panel: fraction of stars that were ejected from the cluster. Bottom panel: Lagrangian radius and the vertical green line marks the moment of core collapse. The time is presented in units of the crossing time of the cluster.

For the case when the mass of the background potential is half the mass in stars, M_ext = 0.5× M_stars the cluster evolves toward core collapse, which now occurs at 1820 t_cross = 16 t_rh (calculated with Eq. (A.14)). The cluster again shows less expansion than the cluster without the external potential (see Fig.B.1) and only one single star is ejected.

For the case when the mass of the external potential is the same as the total mass in stars M_ext = M_stars, the cluster also evolves toward core collapse. However, this is now delayed now until 6553 t_cross as shown by the green vertical line in Fig. B.4.

Fig. B.4. Evolution of a cluster with N = 1000 stars, total mass Mstars = 104 M⊙, Rv = 1.0 pc in a background potential with mass Mext = 104 M⊙. Top panel: fraction of stars that were ejected from the cluster. Bottom panel: Lagrangian radius and the vertical green line marks the moment of core collapse. The time is presented in units of the crossing time of the cluster.

When the number of stars is N = 10 000, the behavior is essentially the same when we include a background potential; the global evolution is delayed due to the increased velocity of the stars. First, we present the case when the mass of the external potential is low M_ext = 0.1 × M_stars compared to the total mass in stars. The cluster also evolves toward core collapse, which now occurs at 3377 t_cross as indicated with a green vertical line in the bottom panel of Fig. B.5, which is ∼14 t_rh (we used Eq. (A.14) to calculate the half mass relaxation time). The mean density inside the 10% Lagrangian radius at this moment is of 4.4 × 10⁵ M_⊙ pc⁻³, which is higher than the mean density at core collapse for the cluster without a background potential. Moreover the cluster in general shows less expansion when including the external potential.

Fig. B.5. Evolution of a cluster with N = 10 000 stars, total mass Mstars = 104 M⊙ and Rv = 1.0 pc under the influence of a background potential with Mext = 103Mext. Top panel: fraction of stars that were ejected from the cluster, which is negligible since only one star escaped the cluster. Bottom panel: Lagrangian radius and the vertical green line marks the moment of core collapse. The time is presented in units of the crossing time of the cluster.

When the mass of the external potential is half the total mass in stars, that is, M_ext = 0.5 × M_stars, the core collapse is delayed until 13 519 t_cross as marked by the green vertical line in Fig. B.6 which is 16 t_rh (we used Eq. (A.14) again to calculate the half mass relaxation time) and during this moment the mean density inside the 10% Lagrangian radius is of 5.5 × 10⁵ M_⊙ pc⁻³. Our simulation is not long enough to see the expansion of the outer parts. however we do expect even less expansion than for the cluster with M_ext = 0.1 × M_stars and again until this time only 1 star has been ejected from the cluster, this was also found in Leigh et al. (2013a); Leigh et al. (2014).

Fig. B.6. Evolution of a cluster with N = 10 000 stars, total mass Mstars = 104 M⊙ and Rv = 1.0 pc under the influence of a background potential with Mext = 5 × 103Mext. Top panel: fraction of stars that were ejected from the cluster, which is negligible since only one star escaped the cluster. Bottom panel: Lagrangian radius and the vertical green line marks the moment of core collapse. The time is presented in units of the crossing time of the cluster.

All Tables

All Figures

	Fig. 1. Evolution of a cluster with N = 10 000 stars with a total stellar mass of Mstars = 104 M⊙ and virial radius Rv = 1.0 pc. Top panel: fraction of ejections from the cluster and bottom panel: Lagrangian radii of the cluster. The vertical green line marks the moment of core collapse.
In the text

	Fig. 2. Evolution of a cluster with N = 10 000 stars with a total stellar mass of Mstars = 104 M⊙ and a virial radius of Rv = 1.0 pc in the center of an external analytic potential with Mext = 104 M⊙. Top panel: fraction of ejections from the cluster and bottom panel: Lagrangian radii of the cluster. The vertical green line marks the moment of core collapse.
In the text

Fig. 3. Core collapse time tcc divided by the crossing time tcross for clusters with N = 1000 (blue circles) and N = 10 000 (red circles) as a function of q = Mext/Mstars. We also plotted two relations that scale as (1 + q)4 which fit our data and therefore suggest that if the core collapse time is proportional to the relaxation time of the cluster, then the relaxation time scales as (1 + q)4, which is the same relation that we derive in Appendix A.

In the text

	Fig. 4. Number of binary systems as a function of time and for different values for the mass of the external potential Mext for simulations with N = 1000 and up to 8000 crossing times. The vertical lines mark the moments of core collapse and coincides with the onset of the formation of binary systems.
In the text

Fig. 5. Distribution of the semi-major axis of the binary population for the set of simulations listed in Table 1. The data clearly shows that increasing the mass of the external potential the formed binaries tend to be more tightly bound. Despite the fact that a low mass of the external potential (Mext <  1.0 ×  Mstars) may promote the formation of wider binaries, these systems tend to be short-lived due to the higher velocities of all the stars. This figure is only intended to show the distribution of the semi-major axis and not the total number of binaries that formed.

In the text

Fig. 6. Evolution of a cluster with N = 1000 stars, total mass Mstars = 104 M⊙, Rv = 0.14 pc, and Rstar = 100 R⊙. Uppermost panel: number of mergers as a function of time, the fraction of binary, hyperbolic, and chain mergers in the second panel along with the fraction of stars that escape from the cluster. Third panel: 10, 50 and 90% Lagrangian radii, and bottom panel: mass growth of the most massive object that formed through mergers in the cluster.

In the text

Fig. 7. Evolution of a cluster with N = 1000 stars, stellar mass Mstars = 104 M⊙, Rv = 0.14 pc, and Rstar = 100 R⊙ in the center of an external potential with Mext = 104 M⊙. Uppermost panel: number of mergers as a function of time, the fraction of binary, hyperbolic, and chain mergers in the second panel, along with the fraction of stars that escape from the cluster. Third panel: 10, 50, and 90% Lagrangian radii, and bottom panel: mass growth of the most massive object that formed through mergers in the cluster. The drop in the 50% Lagrangian radius toward the end of the simulation suggests that the central object may grow up to 5000 M⊙.

In the text

Fig. 8. Example of a Gaussian fit to the number of mergers with the central star for simulation number 45 listed in Table 2. Top panel: number of mergers with the most massive star as a function of time (black solid line) along with the Gaussian function that we used to estimate the total number of mergers at different times. The parameters tdelay and tduration are also shown in the figure. Bottom panel: mass growth of the central star whose rapid growth coincides with the peak in the number of mergers.

In the text

Fig. 9. Mean mass of the most massive star normalized by its initial mass, as a function of N and Rstar for clusters without an external potential evolved until 2000 tcross = 15.6 Myr. The symbols represent the mean values obtained from six simulations and the solid lines show the results of the calculations with our model presented in Eq. (8) using t = 15.6 Myr and the parameters obtained from Eqs. (9)–(12).

In the text

Fig. 10. Mean mass of the most massive star normalized by its initial mass, as a function of N and Rstar for clusters without an external potential evolved until 3000 tcross = 16.5 Myr, except for the point at N = 5000 and Rstar = 20 R⊙, which are simulations that were evolved until 4000 tcross = 22 Myr. The symbols represent the mean values obtained from six simulations and the solid lines show the results of the calculations with our model presented in Eq. (8) using t = 16.5 Myr (t = 22 Myr for N = 5000 and Rstar <  50 R⊙) and the parameters obtained from Eqs. (13)–(16).

In the text

Fig. 11. Parameters of the Gaussian fit to the number of mergers with the central star for simulations of clusters without an external potential. The data points were obtained after the fit to the combined data of six simulations. Upper panel: values of the normalization factor A for different values of N and Rstar along with our model from Eq. (9). These data points show that the number of mergers increases with both N and Rstar as expected. Bottom panel: different values of tdelay in addition to the different values of N and Rstar along with the fitting function presented in Eq. (10). These points show that tdelay decreases with Rstar, but it increases with N.

In the text

	Fig. 12. tduration as a function of the number of stars N for clusters without an external potential. We found no clear correlations between tduration and Rstar. The solid blue line shows the fitting function from Eq. (11).
In the text

Fig. 13. Parameters of the Gaussian fit to the number of mergers with the central star for simulations of clusters with an external potential. The data points were obtained after the fit to the combined data of six simulations. Upper panel: values of the normalization factor Aext for different values of N and Rstar along with our model from Eq. (8). These data points show that the number of mergers increases with both N and Rstar as expected. Middle panel: different values of tdelay, ext in addition to the different values of N and Rstar along with the fitting function presented in Eq. (14). These points show that tdelay, ext decreases with Rstar, but it increases with N. Additionally, these values are larger than the values of tdelay for clusters without an external potential (see Fig. 11). Lower panel: values of tduration, ext along with the fitting function from Eq. (15). We clearly see a positive correlation with N and a negative correlation with Rstar.

In the text

Fig. 14. Fraction of stars fmerger that merge with the central star in a cluster with a total stellar mass of Mstars = 104 M⊙ and a virial radius of Rv = 0.14 pc after 1 Myr (left panel), 2 Myr (middle panel), and 10 Myr (right panel) for different combinations of the number of stars N and the stellar radius Rstar calculated using our model presented in Eq. (8). We can see that when the collisions are maintained for short periods of time, that is, 1 Myr, a massive object with ∼fmergerMstars = 0.15 × 104 M⊙ = 1500 M⊙ is formed in the cluster center only when N is small (given that tdelay becomes larger than the integration time for large N) or large Rstar. If the time limit is long, that is, 10 Myr, most of the clusters will form stars with ∼0.3 Mstars = 3000 M⊙.

In the text

Fig. 15. Fraction of stars fmerger that merge with the central star in a cluster with a total stellar mass of Mstars = 104 M⊙, a virial radius of Rv = 0.14 pc, and an external potential with a mass of Mext = Mstars after 1 Myr (left panel), 2 Myr (middle panel), and 10 Myr (right panel) for different combinations of the number of stars N and the stellar radius Rstar calculated using our model presented in Eq. (8). We can see that when the collisions are maintained for short periods of time, that is, 1 Myr, an object that is not so massive with ∼fmergerMstars = 0.1 × 104 M⊙ = 1000 M⊙ is formed in the cluster center only when N is very small (given that tdelay becomes larger than the integration time for large N) or when Rstar is large. If the time limit is 2 Myr (middle panel), we expect that clusters with a small N experience core collapse and form a massive central star. However, if the number of stars is large and Rstar is also large, given the higher velocity of the stars in those clusters compared to clusters without the external potential, there is a larger probability for hyperbolic mergers to occur before core collapse, which we see on the right-hand side of middle panel. If the time limit is long, that is, 10 Myr, most of the clusters form stars with ∼0.3−0.4 Mstars = 3000−4000 M⊙ and even 0.7 Mstars = 7000 M⊙ if N is large enough.

In the text

	Fig. A.1. Trajectory of a star with mass M that passes close to another star of mass m, which causes a small deflection on the velocity vector V.
In the text

	Fig. B.1. Evolution of a cluster with N = 1000 stars, total mass Mstars = 104 M⊙ and Rv = 1.0 pc. Top panel: fraction of stars that were ejected from the cluster. Bottom panel: Lagrangian radius and the vertical green line marks the moment of core collapse. The time is presented in units of the crossing time of the cluster.
In the text

	Fig. B.2. Evolution of a cluster with N = 1000 stars, total mass Mstars = 104 M⊙ and Rv = 1.0 pc in a background potential with mass Mext = 103 M⊙. Top panel: fraction of stars that were ejected from the cluster. Bottom panel: Lagrangian radius and the vertical green line marks the moment of core collapse. The time is presented in units of the crossing time of the cluster.
In the text

	Fig. B.3. Evolution of a cluster with N = 1000 stars, total mass Mstars = 104 M⊙, Rv = 1.0 pc in a background potential with mass Mext = 5 × 103 M⊙. Top panel: fraction of stars that were ejected from the cluster. Bottom panel: Lagrangian radius and the vertical green line marks the moment of core collapse. The time is presented in units of the crossing time of the cluster.
In the text

	Fig. B.4. Evolution of a cluster with N = 1000 stars, total mass Mstars = 104 M⊙, Rv = 1.0 pc in a background potential with mass Mext = 104 M⊙. Top panel: fraction of stars that were ejected from the cluster. Bottom panel: Lagrangian radius and the vertical green line marks the moment of core collapse. The time is presented in units of the crossing time of the cluster.
In the text

Fig. B.5. Evolution of a cluster with N = 10 000 stars, total mass Mstars = 104 M⊙ and Rv = 1.0 pc under the influence of a background potential with Mext = 103Mext. Top panel: fraction of stars that were ejected from the cluster, which is negligible since only one star escaped the cluster. Bottom panel: Lagrangian radius and the vertical green line marks the moment of core collapse. The time is presented in units of the crossing time of the cluster.

In the text

Fig. B.6. Evolution of a cluster with N = 10 000 stars, total mass Mstars = 104 M⊙ and Rv = 1.0 pc under the influence of a background potential with Mext = 5 × 103Mext. Top panel: fraction of stars that were ejected from the cluster, which is negligible since only one star escaped the cluster. Bottom panel: Lagrangian radius and the vertical green line marks the moment of core collapse. The time is presented in units of the crossing time of the cluster.

In the text