Issue 
A&A
Volume 620, December 2018



Article Number  A70  
Number of page(s)  12  
Section  Galactic structure, stellar clusters and populations  
DOI  https://doi.org/10.1051/00046361/201833854  
Published online  03 December 2018 
The hunt for selfsimilar core collapse
^{1}
Astronomical Institute of Charles University, Prague, Czech Republic
email: pavlik@sirrah.troja.mff.cuni.cz.
^{2}
Observatory and Planetarium of Prague, Prague, Czech Republic
Received:
13
July
2018
Accepted:
13
August
2018
Context. Core collapse is a prominent evolutionary stage of selfgravitating systems. In an idealised collisionless approximation, the region around the cluster core evolves in a selfsimilar way prior to the core collapse. Thus, its radial density profile outside the core can be described by a power law, ρ ∝ r^{−α}.
Aims. We aim to find the characteristics of core collapse in Nbody models. In such systems, a complete collapse is prevented by transferring the binding energy of the cluster to binary stars. The contraction is, therefore, more difficult to identify.
Methods. We developed a method that identifies the core collapse in Nbody models of star clusters based on the assumption of their homologous evolution.
Results. We analysed different models (equal and multimass), most of which exhibit patterns of homologous evolution, yet with significantly different values of α : the equalmass models have α ≈ 2.3, which agrees with theoretical expectations, the multimass models have α ≈ 1.5 (yet with larger uncertainty). Furthermore, most models usually show sequences of separated homologous collapses with similar properties. Finally, we investigated a correlation between the time of core collapse and the time of formation of the first hard binary star. The binding energy of such a binary usually depends on the depth of the collapse in which it forms, for example from 100 kT to 10^{4} kT in the smallest equalmass to the largest multimass model, respectively. However, not all major hardenings of binaries happened during the core collapse. In the multimass models, we see large transfers of binding energy of ∼10^{4} kT to binaries that occur on the crossing timescale and outside of the periods of the homologous collapses.
Key words: methods: numerical / galaxies: clusters: general / binaries: general
© ESO 2018
1. Introduction
Twobody relaxation states that any star moving through a field of stars is decelerated by the force of dynamical friction (Chandrasekhar 1943), which is proportional to the mass of the moving star and inversely proportional to its velocity squared. This force is also responsible for mass segregation in multimass systems. Since the heat capacity of a selfgravitating system in virial equilibrium is negative (e.g. LyndenBell & Wood 1968; Binney & Tremaine 1994), the central region of a star cluster should contract over time. Consequently, the central parts relax more quickly than the halo and the velocity distribution in the centre is almost Maxwellian (Larson 1970). Within the thermodynamic framework, the core is supposed to collapse, reaching infinite density and kinetic temperature in a finite time (also known as the gravothermal catastrophe). Unlike continuum models, in Nbody models (and real star clusters) this sequence is prevented by the presence of existing or newly formed binary stars in the core, whose ability to efficiently expel other stars via threebody interactions cools the core (e.g. Aarseth 1972; Hut 1983; Fujii & Portegies Zwart 2014; O’Leary et al. 2014). The cluster core gradually shrinks towards the collapse and then expands rapidly (so called core bounce). Thus, the event of core collapse may be indirectly observed but its exact time is no longer well defined.
LyndenBell & Eggleton (1980) showed that the evolution of a spherically symmetric collisionless system prior to core collapse should be selfsimilar (homologous), that is its density evolves with respect to the radius and time according to the scaling relation
Here, ρ_{c} is the core density, ρ_{⋆} is a dimensionless structure function, and the radius is described using an enclosed mass m and a scaling factor r_{⋆},
where r_{c} stands for the core radius and is the core mass. The homologous solution implies that the internal structural scaling has an exponent α that remains temporally invariant. As it must also satisfy smoothness conditions for ρ_{⋆}(r_{⋆}) and normalisation, generally α = const. (LyndenBell & Eggleton 1980, also e.g. Penston 1969). The core radius then depends on the core density as and the temporal evolution of the core radius before the time of core collapse, t_{cc}, is
According to LyndenBell & Eggleton (1980), the radial density profile may be approximated by a doublebroken powerlaw function with the logarithmic density gradient defined as
which is equal to zero in the cluster core and reaches α asymptotically. In an intermediate region above r_{c}, the logarithmic density gradient has to be larger than α to compensate for the missing mass in the core (see the slopes a_{I}, a_{II}, and α in Fig. 1). LyndenBell & Eggleton (1980) found α ≈ 2.208. Further works based on either isotropic (Cohn 1980) or anisotropic models (Takahashi 1995) in a Fokker–Planck approximation led to a slightly different value, α ≈ 2.23.
Fig. 1. Schematic plot of the radial density profile in an Nbody cluster during the core collapse with four different values of the logarithmic density gradient, a_{I − IV}. The dotted line has a slope equal to α and shows the asymptotic solution of LyndenBell & Eggleton (1980). The break radius r_{I} is identified with the core radius, r_{III} roughly corresponds to the halfmass radius, and r_{k} is the cluster radius. The notation is the same as in Eq. (15) that we used for the fitting of our numerical models. 

Open with DEXTER 
In Nbody star clusters, the selfsimilar solution is not infinite. In this case, distant parts of the halo tend to a Maxwellian distribution of velocities on the relaxation timescale. This changes the logarithmic density gradient in the halo to a = 3.5 (Spitzer & Hart 1971). Hence, we expect the cluster’s radial density profile to be approximated by a triplebroken power law as indicated by a solid line in Fig. 1.
Radial density profiles of numerical Nbody models (e.g. Giersz & Heggie 1994; Makino 1996b; Baumgardt et al. 2003) show a good agreement with the expectations above. In this paper we go beyond those works in several aspects: (i) computational resources available nowadays enable us to integrate and analyse hundreds of realisations of models consisting of tens of thousands of particles; (ii) besides equalmass models, we evaluate models with a Salpeter mass function and; (iii) we not only analyse the density profile of the cluster at the time of core collapse but also test the hypothesis of a homologous evolution in time which, among other things, allows us to formulate a new method for determining the time of core collapse.
2. Numerical models
The numerical models that we are working with may be scaled arbitrarily, so we express all numbers in this paper in Hénon units (unless we specify otherwise) that are used in Nbody integrations, that is, G = M_{tot} = R_{virial} = 1. The total energy is and the crossing time is in those units.
We studied several Nbody models of star clusters: two equalmass (e10k, e50k) and three multimass (m20k, m100k, n100k) with a Salpeter (1955) slope of −2.35 for the initial mass function in the mass range listed in Table 1. We generated initial conditions of all models using the plumix code (Šubr et al. 2008; Šubr 2012) without initial mass segregation, that is with a Plummer (1911) distribution function. The dynamical evolution was followed using the collisional nbody6 code (Aarseth 2003).
Parameters and notation of the numerical models analysed in this paper.
Being concerned only with dynamical effects, our star clusters were modelled as isolated systems, and we did not consider any external potential or tidal forces from the Galaxy. Stars were represented by point masses and we did not include any primordial binaries or stellar evolution. Collisions of stars were disabled and escaping stars were not removed from any simulation. All models are pure Nbody and neither gas nor dust components were included. In other words, this work is an exploration of the dynamics of isolated, ideal selfgravitating systems.
3. Methods
3.1. Time of core collapse
One definition of the time of core collapse is the moment of the minimum core radius. However, when the cluster core shrinks, its binding energy rises and all stars that are moving inward or outward have increasingly greater speeds. The core becomes unstable and subject to random pulsations (Makino 1996a), and the definition of the core is not obvious. Consequently, the time of its minimum radius is not clearly defined. An alternative method, which we will discuss later, is related to the formation of the first hard binary system (Fujii & Portegies Zwart 2014). In this section, we present a third method that is based on analytic predictions of the selfsimilar evolution of a star cluster before the core collapse.
The evolution of star clusters is often described by means of the Lagrangian radii (e.g. see Figs. A.1–A.4), which are defined as the radii of concentric spheres containing fixed mass fractions. This definition is obviously dependent on the choice of the origin of the coordinate system because at each time step, stars that occupy these spherical regions are reordered by their current radial distance from the origin. In our models, we calculate the Lagrangian radii from the cluster’s density centre provided by nbody6 (i.e. a numerically improved method from Casertano & Hut 1985, which is a well established and commonly used approach).
In a homologous solution, the evolution of the Lagrangian radii is given by Eq. (2). At the moment of their minima, the Lagrangian radii are related to r_{c} by a single constant of proportionality. The minimum of every Lagrangian radius curve is given by
After substituting r in (5) from Eq. (2), we get
where a dot represents and, using homology, m_{⋆} = m/m_{c}(t). Following the same argument as LyndenBell & Eggleton (1980) that each side of this equation depends on a different parameter (t and m_{⋆}, respectively), both sides must be constant with respect to these parameters. Hence, the scaling functions m_{⋆} and r_{⋆} have the same value for all minima of the Lagrangian radii. This implies that the radius representing the minimum (hereafter r_{m}) is a fixed multiple of r_{c} at the time of its minimum (denoted as t_{m}). Together with Eq. (3), we have
Consequently, parameters of the homologous solution, t_{cc} and α, are obtainable from fitting function (7) to the sequence of minima of the Lagrangian radii in numerical models.
3.1.1. Data processing and fitting
As mentioned above, the cluster core pulsates with unpredictable frequency and amplitude while collapsing. Fluctuations of the unstable core also propagate, by definition, in the Lagrangian radii. They are most visible for low mass fractions, that is, inner Lagrangian radii. These fluctuations arise from random motions of stars in the cluster core. We are, however, interested in general trends that are driven by the cluster thermodynamics. In order to reveal them and, in particular, to find the corresponding minima of the Lagrangian radii for Eq. (7), we reduced the fluctuations by smoothing each Lagrangian radius timewise. We tried a moving average and an algorithm from Savitzky & Golay (1964), which is a method with a higher order polynomial. The latter proved better in reducing the fluctuations while preserving the global evolution of the Lagrangian radii, so we will focus only on that one.
We found the minima of n smoothed Lagrangian radii that provide us with a set of points
Assuming that these points of minima satisfy the homologous solution described by Eq. (7), we may find t_{cc} by fitting them with function
via the parameters t_{cc}, α, and a proportionality constant A. This function may only be used for t_{m} ≤ t_{cc}, which makes t_{cc} not only a fitting parameter but also an upper bound for t_{m}. It is evident that, for example Eq. (9) cannot be used with a fitting method that requires fixed limits; the solution in such cases is to introduce a Heaviside function
and then redefine Eq. (9) to
where t_{m} is formally no longer restricted by the choice of bounds.
The best fit is determined by the lowest value of a sum of the squared absolute deviations in radii, which we call , defined as
As a minimising method, we applied a genetic algorithm by Storn & Price (1997) ^{1}.
3.1.2. Limits
The Storn–Price minimisation algorithm requires fixed lower and upper limits for each parameter to generate an initial random sample of seeds in the parameter space. The range of possible values of α must satisfy the condition α ≠ 6 that comes from the exponent in Eq. (3). According to the analytic solution, α is supposed to reach a value between 2 and 2.5 (cf. LyndenBell & Eggleton 1980), hence we set the upper limit to α < 6. Assuming that the density is a decreasing function of distance from the centre, we set the lower limit to α ≥ 0.
The points of minima (8) are supposed to represent a decreasing powerlaw function, so we do not expect the time of core collapse to happen sooner than the earliest minimum of the Lagrangian radii we fitted. Therefore, the lower limit of t_{cc} was set to min(t_{m, i}). The upper limit is a value that is reasonably large, so it will not affect the results, that is, 2max(t_{m, i}) − min(t_{m, i}). The range of the proportionality constant A was set wide enough not to restrict the fit, for example from −2 to 2, whereas the typical values of A were around 0.008 in model e10k and around 0.04 in model m20k.
3.2. Density profile at core collapse
To evaluate the radial density profile of a star cluster during the core collapse, we linearly interpolated each Lagrangian radius to the determined t_{cc} from the nearest values around that time. Hereafter, we denote the ith Lagrangian radius at core collapse by r_{I}. So, at t_{cc}, we have a set of Lagrangian radii
where r_{1} corresponds to the lowest mass fraction, r_{k} is the cluster radius (defined here by the Lagrangian radius of 99% of mass), and r_{0} ≡ 0. The mean density contained in the spherical shell between every pair of consecutive Lagrangian radii (assuming spherical symmetry) is given by
where m_{I} is the mass percentage (i.e. the mass in Hénon units) contained in a sphere of radius r_{I}. By definition, m_{0} = m(r_{0})≡0. Finally, we have a set of points {(r_{I},ρ_{I})} representing the radial density profile of a star cluster at the time of core collapse.
3.2.1. Fitting
We fit a triplebroken powerlaw function (see Fig. 1)
to the set {(r_{I},ρ_{I})}. The variables r_{I}, r_{II}, and r_{III} denote the radii where this powerlaw function breaks, and r_{k} is the cluster radius (in our case the Lagrangian radius of 99% of mass). Following the discussion in Sect. 1, during the core collapse we may also use the fit to determine the core radius r_{c} ≡ r_{I}. Requiring function (15) to be continuous, the proportionality constants b_{II − IV} must satisfy the conditions
and
where we defined b ≡ b_{I} as it is now the only proportionality constant to fit. We again applied the Storn–Price algorithm to minimise the sum of the squared absolute deviations in the radial density
using parameters a_{I − IV}, r_{I − III} and b.
3.2.2. Limits
The radial density profile (15) is expected to be a strictly decreasing function above the core radius. In addition, the slope a_{III} is supposed to have a value close to α. Therefore, we defined common boundaries of a_{II − IV} ∈ [0, 6]. According to LyndenBell & Eggleton (1980), the radial density profile is flat in the core (a_{I} ≈ 0), but working with discrete data, binning has a great impact on the calculation of the Lagrangian shells. Therefore, we allowed for the density also having a positive slope in the core, giving it a generous range a_{I} ∈ ( − 6, 6). We took the inner and outermost Lagrangian radii for the limits of r_{I − III}, so r_{I − III} ∈ (r_{1}, r_{k}). We also used an implicit condition that r_{I} < r_{II} < r_{III}. A typical value of the proportionality constant was for example b ≈ 1.6 in model e10k and b ≈ 0.8 in model m20k. Hence, we set the range of b in all four models reasonably wide (e.g. from −5 to 5) not to restrict the fits.
4. Results
The method we described above was applied to all of our models (e10k, e50k, m20k, m100k, and n100k). First, we discuss the noise reduction. Then, for each model individually, we focus on finding the minima, and fitting the time of core collapse and radial density profile. Finally, we compare our results with another method for estimating the time of core collapse via the formation of hard binary stars.
4.1. Data preparation
As we mentioned above in Sect. 3.1, oscillations of the unstable collapsing core make it impossible to seek systematic evolution of the Lagrangian radii using the raw data. We reduced those fluctuations using the (Savitzky–Golay) algorithm with a second order polynomial, and the window width of the order of ten crossing times for each model. Results of this procedure are shown in the upper panels of Figs. A.1–A.4, where we plot the evolution of one arbitrarily chosen realisation of models e10k, e50k, m20k, and m100k. The effect of this data filtering is most visible in the inner Lagrangian radii of both plotted multimass models, where the noise is very high due to a small number of massive stars dominating the cluster’s central region.
Alternatively, averaging the Lagrangian radii over many realisations (bottom four panels of Figs. A.1–A.4) also substantially reduces random fluctuations. For the multimass model, it renders the time of core collapse clearly visible as the global minimum of the inner Lagrangian radii. This may not be so obvious in a single realisation, where the global minima of the inner Lagrangian radii (either smoothed or not) may occur at later times (see the case presented in Figs. A.3 and A.4). However, the advantage of the (Savitzky–Golay) filtering is that it not only allows us to identify t_{cc} in individual realisations, based on the position of the local minima of the inner Lagrangian radii, but it also reveals longterm oscillations of the core region, which could be related to gravothermal oscillations in more populous models (Makino 1996b).
4.2. Equalmass models
In the case of model e10k, we have integrated each realisation twice, first for an overview of the global evolution, as shown in Fig. A.1. When the contraction of the inner cluster region started, we reintegrated the models once more with a more frequent output until the contraction ended in order to have high resolution data for the following analyses. Filtering these outputs provided us with a clear single global minimum of each Lagrangian radius. The set (8) was then constructed from these minima of the Lagrangian radii of mass fractions m_{i, i ≥ 1} ∈ [0.001, 0.03] with a step of 0.001: those points are plotted for an arbitrarily chosen realisation in the top left panel of Fig. A.5. In this particular realisation, the fitting procedure described in Sect. 3.1 gives t_{cc} ≈ 2317 and α ≈ 2.215. After the evaluation of all one hundred realisations, we found that the dispersion of the fitted t_{cc} indicates that the core collapse happened at a comparable time 2297 ± 52 (see Table A.2). This result is consistent with Fujii & Portegies Zwart (2014) who argued that the time of core collapse should be characteristic for any given model, depending only on the ratio between the mass of the most massive star and the mean stellar mass. Our fit of the powerlaw index α = 2.33 ± 0.02 (see Table A.2) also agrees with theoretical expectations (LyndenBell & Eggleton 1980; Cohn 1980; Takahashi 1995).
At the time of core collapse, we constructed the radial density profile using the method described in Sect. 3.2. The density profile (see the bottom left panel of Fig. A.5 for an example) is in accord with our expectations: it is almost flat in the core, steeper around the core (slightly shallower at larger radii), and the steepest in the halo (see Table A.3 for the slopes and Table A.4 for the radii where the density profile breaks). The fitted slope a_{III} ≈ 2.321 from (15) agrees with previous Nbody simulations (e.g. Baumgardt et al. 2003). The value of a_{III} is nearly identical to the powerlaw index α ≈ 2.33 obtained by fitting the Lagrangian radii. The near equality of α and a_{III} indicates that the equalmass model really evolved selfsimilarly before the core collapse. The outermost logarithmic density gradient, a_{IV} ≈ 3.4, is also in a good agreement with the prediction of Spitzer & Hart (1971).
After the first major contraction of the core region in model e10k, the inner Lagrangian radii show at least one postcollapse oscillation before the end of integration (a similar behaviour was described also e.g. by Makino 1996b). In such a small cluster, which barely exceeds the number of particles needed to form an unstable core (cf. Goodman 1987), post collapse oscillations are not very deep and thus hard to detect. With a higher number of stars in a model, the core density increases and the collapse will be deeper (e.g. Hut 1996; Fujii & Portegies Zwart 2014). Due to these reasons, we have integrated an additional model (e50k) of a more massive equalmass star cluster, where those gravothermal oscillations are more easy to quantify. In order to determine the time of core collapse and the time of the subsequent contractions, we used the first and second time derivatives of each Lagrangian radius (calculated with finite differences) to identify the points of local minima. To eliminate false detection, we took a moving window with the same width as in the smoothing (of the order of ten crossing times) and kept only the deepest minimum in that window. Then we treated each of these sequences separately. For the region around the core, we used mass fractions from 0.001 to 0.020 with a step of 0.001 (see the first sequence in the top right panel of Fig. A.5; the plotted realisation is one of the few that had two postcollapse oscillations in the integration window). At each minimum, we constructed the radial density profile of the cluster.
The first major contraction in e50k occurred at time t_{cc} ≈ 9350 (see Fig. A.2); the subsequent contractions appear at a slightly different time and their depth varies in individual realisations. The homologous index reads α ≈ 2.24 at the core collapse in all realisations. This value is in agreement with LyndenBell & Eggleton (1980) and comparable with the slope a_{III} (see Table A.3). The radial density profile is in good agreement with model e10k and the predicted slopes (cf. LyndenBell & Eggleton 1980; Spitzer & Hart 1971). We also fitted the time of the first postcollapse oscillation in e50k, which was present in all realisations of this model (e.g. see the second sequence in the top right panel of Fig. A.5), and its corresponding radius density profile (the green curve in the bottom right panel of Fig. A.5). The plots and the results printed in Tables A.2 and A.3 show that in this case α ≈ 2.34 and a_{III} ≈ 2.37. Thus, we claim that the first postcollapse oscillation is homologous, as well as the core collapse.
4.3. Multimass models
The inner Lagrangian radii of both m20k and m100k models show multiple lowfrequency waves that represent several contractions of the core region in the time span of tens of crossing times. Those contractions are unevenly distributed in time in the individual realisations, that is they are not visible in the averaged radii (cf. upper and lower panels in Figs. A.3 and A.4). Only the first contraction, which corresponds to the core collapse, seems to happen at the same time in all realisations (t_{cc} ≈ 53 in m20k and t_{cc} ≈ 116 in m100k). Because the depths of these contractions vary, we used derivatives of the Lagrangian radii to determine the sequences of minima, as in the case of e50k. In 80% of realisations we found two and in 45% of realisations three clear consecutive minima. Due to a greater effect of binning, in this model we used larger mass fractions m_{i, i ≥ 1} ∈ [0.005, 0.03] with a step of 0.005 (see the top panels of Fig. A.6) when constructing each set (8).
We fitted each of these sequences separately (see an example of this procedure in the right panels of Fig. A.6; the corresponding mean values of t_{cc} and α are listed in Table A.2). Although the points of minima are not as ordered as in the equalmass models, the first contraction is well defined across all realisations and its deviation is rather small. The mean times of the second and third minima have higher uncertainties and their 1σ intervals overlap. This also proves that those minima are unevenly distributed in time across all realisations. The derived values of α, about 1.5, are significantly different from the theoretical predictions (α ≈ 2.2) but within their uncertainties, they are the same in all three minima.
For each estimated time of core collapse, we constructed the radial density profile (see Table A.3). Qualitatively, it follows our expectations: almost zero in the centre, then a steep slope followed by a shallower one, and the steepest in the halo. However, the values of all slopes are different from those in the equalmass models. Greater uncertainties that we see in the fits are inevitable due to a small number of particles that produce the Lagrangian radii of small mass fractions; each radius is highly influenced by any massive star that passes through the central region. Further, we note that the halo has a much steeper profile than predicted analytically. Nevertheless, the value of a_{III} ∈ (1.6, 1.7) (see Table A.3) is similar in all minima and compatible with α within its uncertainties. This result indicates that the multimass clusters could evolve selfsimilarly too, albeit with different parameters than equalmass or analytic models.
The detected postcollapse contractions of the cluster core most probably do not represent gravothermal oscillations because they are too separated in time and the systems are too small to form an unstable core (e.g. Breen & Heggie 2012a,b). We rather refer to them as homologous collapses; they look almost the same and have comparable (and perhaps selfsimilar) properties. This leads us to conclude that in a more complicated system (e.g. a real star cluster, which we are well aware our models are not) simply analysing its postcollapse dynamical structure may not be enough to distinguish which collapse it has already sustained.
The depth of core collapse and its subsequent oscillations depends on the number of massive stars in the model and the ratio between the mass of the most massive star and the total mass of the cluster (Breen & Heggie 2012a,b). In both m20k and m100k, this ratio was of the order of 10^{−3}. In order to approach the depth of core collapse of our equalmass models (where this ratio is 10^{−4} or 10^{−5}), we made an additional multimass model (n100k) with the same initial mass function slope as models m20k and m100k, but with a slightly modified range of masses to acquire a ratio of 10^{−4} (see Table 1). Even in this system, the collapse was not as deep as we expected. Fluctuations of the inner Lagrangian radii made our method for finding selfsimilar core collapse inefficient as can be seen for instance in the top panel of Fig. A.7 where our method detected two equivalent collapses of which only one is a potential core collapse. Based on the radial profile constructed at the times of minima, we got the powerlaw slope a_{III} ≈ 1.9 but the index α was different in most cases, ranging from 1.0 to almost 3.2, and in some cases, the sequence of minima was increasing instead of decreasing.
Selfsimilar evolution is a feature of collisionless systems. Although we found traces of it in some collisional systems, in others where for example close encounters are more dominant, we did not. Based on our results, we are unable to make a general statement on whether multimass models evolve selfsimilarly prior to the core collapse or not. In order to evaluate selfsimilar evolution in star clusters with a mass function, larger models with higher N or a modified approach would be needed.
4.4. Binary star formation
The binding energy of a binary star (composed of stars with masses m_{a} and m_{b}, with a semimajor axis d) is
We express the binding strength of a binary star in terms of kT, where k is the Boltzmann constant and T is the kinetic temperature, which relates to Nbody variables (in Hénon units) by
Once a binary star becomes sufficiently tightly bound (E_{bin} exceeds several kT) it has a very low probability of being destroyed due to fewbody interactions (Heggie 1975; Goodman & Hut 1993). To acquire that amount of energy, the pair must live in a dense environment, such as a collapsing cluster core (Tanikawa et al. 2012).
It has been suggested (e.g. Fujii & Portegies Zwart 2014) that the time, t_{bin}, of the first appearance of a hard binary star with a critical energy
where m_{max} and ⟨m⟩ are the maximum and mean mass, respectively, identifies the core collapse. In the following, we test this hypothesis using our independent method for finding the time of core collapse.
For our model e10k, where m_{max} = ⟨m⟩, Eq. (22) gives E_{lim} = 10 kT. For that value of energy, we found a very good correlation between t_{bin} and t_{cc} (coefficient 0.884, see Table A.1 and left panel of Fig. A.10). A slightly better correlation coefficient of 0.909 was found for E_{lim} = 100 kT, which indicates that the process of energy transfer into binaries is very quick, as was also pointed out by Tanikawa et al. (2012). During the process of core collapse (i.e. in the range of a few crossing times around t_{cc}) a typical binary star acquires energy ΔE_{bin} of a few hundreds or even a thousand kT (see the histogram in Fig. A.8), which is at least one order of magnitude greater than the binding energy which is supposed to identify the core collapse.
Equation (22) gives E_{lim} ≈ 750 kT for the m20k model. Following the same method as in e10k, we calculated the correlation coefficients for various binding energies of binary stars in the first collapse. Our choice of limiting binding energy comes from the histogram in Fig. A.9, which shows the increase of binding energy of the hardest binary during each collapse. We found equivalent correlations for a range of energies between 750 kT and 1250 kT (compare the results in Table A.1 and Fig. A.12), which also supports the idea brought by Tanikawa et al. (2012).
In the case of model n100k, the limiting energy for a hard binary is E_{lim} ≈ 100 kT. In this model, we have detected several binary stars with E_{bin} > E_{lim} existing at t ≈ t_{cc} (see the lower panel of Fig. A.7). The first major contraction in the plotted realisation was driven by several binaries and the second seems to form a 13 000 kT binary star (which is well above E_{lim}). Due to a high uncertainty in determining the time of core collapse in this model, we were unable to calculate the correlation with the time of formation of the first hard binary star.
Let us also point out that rapid hardening of binaries, or formation of new ones, is not unique to the core collapse. We found several occurrences of such events even though we did not find any sign of collapse in the star cluster (e.g. compare the process of binary evolution in models m20k and m100k in the bottom panels of Fig. A.6). On the left hand side, we see the evolution of two hard binaries (plotted in red and blue) in model m20k. The red binary star formed during the first collapse of the core and hardened continuously. In the subsequent homologous collapse, a shortlived binary star (green) emerged and the existing one (red) acquired almost twice its former binding energy. A detailed analysis showed that the disappearance of the green binary was due to a close interaction with the red one at that time. During the third collapse, another hard binary (blue) started to form and harden. The event of large energy change of a binary that occurred out of sync with any collapse is clearly visible on the red evolution track.
In the bottom right panel of Fig. A.6, we show a similar process in one realisation of model m100k. Although the formation of the first binary (red) began earlier than the core collapse determined from the fitting, the first major increase of its binding energy is linked to the first contraction of the core. A detailed look at the output revealed that a new binary (green), which formed at t ≈ 200, interacted at t ≈ 250 with the red binary star. This caused an immediate hardening of the final binary, which follows in red, and interchange of the components. During the second homologous collapse in the plotted realisation, the red binary hardened and it was expelled from the cluster, and a new binary (blue) started to form. It gradually hardened and is likely to be related to the third homologous collapse. Also in this realisation, we see cases of formation of binaries (the green one) or their rapid hardening (the blue binary after time 450) that are out of sync with the collapses of the core.
The correlation between the time of core collapse and the time of the formation of the first binary star in the more massive equal and multimass models (i.e. e50k and m100k) is based on the data from only ten realisations. Hence, characterising a typical energy of a binary star that was promoted by the core collapse is influenced too much by the statistical noise. We do not draw any conclusions from these particular results. Nevertheless, we may observe a similar trend as in the smaller models, that is that the core collapse is linked with the formation of very hard binary stars. In the case of e50k (see Fig. A.11), the best correlation is achieved for E_{lim} between 100 kT and 1000 kT with a coefficient above 0.98. In the case of model m100k, the correlation coefficients vary from 0.47 to 0.78 for E_{lim} above 1500 or 10^{4} kT, respectively (compare also the plots in Fig. A.13). These energies are again well above the estimate from Eq. (22).
We did not make an attempt to evaluate the correlation of formation (or hardening) of binaries with the second or third homologous collapse as it is virtually impossible to distinguish which collapse is responsible for creating a binary (or vice versa) after the system has already collapsed once.
5. Conclusions
We investigated the properties of core collapse in numerical Nbody models of selfgravitating star clusters. For that purpose, we developed a novel method for the identification of the time of core collapse. The method is based on an assumption proposed for analytic models by LyndenBell & Eggleton (1980) that the evolution of the cluster is selfsimilar.
In the case of equalmass models (e10k and e50k), we found a very good agreement with theoretical expectations. Minima of the Lagrangian radii for small mass fractions are aligned according to a powerlaw relation with the powerlaw index close to α ≈ 2.3. At the time of core collapse, the cluster’s radial density profile in the intermediate region between the core and the halfmass radius is well approximated by a power law ρ ∝ r^{−aIII}, with a_{III} ≈ 2.3. The fact that a_{III} ≈ α indicates that the cluster’s evolution matches the selfsimilar solution of LyndenBell & Eggleton (1980). The density profile in the halo is best fitted by a power law with index a_{IV} ≈ 3.4, which is close to the prediction formulated by Spitzer & Hart (1971) for the evolution of halos of Nbody models.
Further, we analysed Nbody models of star clusters with a Salpeter (1955) mass function (m20k and m100k). Using our method, we identified the times of core collapse and determined the radial density profile of the clusters at that moment. We found that the cluster’s evolution and density profile are qualitatively similar to the previous case, although the powerlaw index α has a significantly different value. Specifically, the bestfit value of α for temporal evolution of the inner Lagrangian radii is 1.5, which is nearly identical to the powerlaw index of the radial density profile beyond the cluster core, a_{III} ≈ 1.6. Thus, we conclude that these models show traces of selfsimilar evolution.
We also studied the evolution of a multimass model (n100k) with the same slope of the mass function as m20k and m100k but a higher ratio between the total mass and the most massive star. In terms of selfsimilar evolution, we expected this model to be a “bridge” between the equal and multimass models that we have already discussed. However, there were big differences in the radial profiles across the realisations, caused by random oscillations of the core region. In most realisations, we were unable to successfully fit the minima of the Lagrangian radii and clearly determine the time of core collapse and its homologous properties.
Our results show that analytic predictions on the selfsimilar evolution are valid in the limit of equalmass Nbody systems but cannot be straightforwardly generalised for multimass (i.e. more realistic) star clusters. A further study from both the analytic and numerical point of view is needed to conclude whether multimass systems with a general mass function do undergo selfsimilar core collapse evolution, perhaps with the homologous index dependent on the mass function properties. Any future studies of this topic would certainly benefit from analysing even more populous clusters.
In the case of m20k and m100k as well as in e50k, we found subsequent phases of coherent evolution of the inner Lagrangian radii even after the core collapse. Evolution toward all those minima have similar characteristics and homologous properties (i.e. depth of the core contraction, powerlaw indices α, and the radial density profiles). Therefore, we conclude that they are observationally indistinguishable from each other. The only prominent difference between the first and subsequent homologous collapses is that the time of the first one is well correlated among different realisations, which corresponds to the argument made by Fujii & Portegies Zwart (2014). Our values of t_{cc} are 52.9 ± 8.1 (m20k) and 116 ± 38 (m100k), while for the second and third collapses, for example in m20k, we have found 101 ± 21 and 120 ± 17, respectively (see also Table A.2). A large deviation of the times of subsequent collapses implies that they are smeared out in the plots of the Lagrangian radii averaged over all realisations of the particular models (see Figs. A.3 and A.4). In the case of m20k, we identified two such homologous collapses (including the first one) in 80% and three in 45% of the realisations within the integration time. All realisations of m100k passed at least three homologous collapses.
Finally, we studied the correlation of the time of core collapse, t_{cc}, determined by our method with the formation of dynamical binaries in the cluster. In the case of e10k, we found the best correlation of t_{cc} with the time when the first binary acquired binding energy higher than 100 kT (correlation coefficient of 0.909), yet only a slightly smaller value (≈0.88) was obtained for the correlation with the first occurrence of a binary with E_{bin} > 10 kT. This indicates that (i) the flow of energy toward the binaries is indeed very fast during the core collapse and (ii) the formation of the first hard binary with relatively poorly constrained binding energy may be used to identify the core collapse. In the multimass model we have the best correlations for binding energies between 750 kT and 1250 kT, where the correlation coefficient is in the range from 0.53 to 0.58. Analytic estimates for the binding energy of binaries formed during the core collapse derived by Fujii & Portegies Zwart (2014) give values of 10 kT and 750 kT for the e10k and m20k models, respectively.
Detailed inspection of our models revealed that a large transfer of binding energy from the cluster to binaries occurs not only during the core collapse but also during the subsequent homologous collapses. On the other hand, tracking the binding energy of binaries in our models (m20k in particular) revealed that episodes of large energy transfer are much more numerous than the homologous collapses. In other words, there are common events of formation of (or hardening of existing) binaries that cannot be identified with any homologous collapse. In some cases, these interactions led to a change of binding energy of the order of 10^{4} kT on a timescale shorter than one crossing time, exceeding by an order of magnitude the energy transfer rate related to the homologous collapses.
The formation of the first hard binary star (in a system without primordial binaries) is well correlated with the phase of core collapse in a star cluster. As there are no other hard binaries present in the system, it is a good indicator of this event. After the system has already collapsed once and has produced at least one hard binary star, neither the formation of a new hard binary nor a large transfer of binding energy into existing ones can be considered as an indicator of the subsequent homologous contractions. From that perspective, it would be intriguing to examine homologous properties and binary evolution during the core collapse in systems containing a primordial binary population.
Acknowledgments
First and foremost we thank Douglas Heggie for introducing us to this topic and for his selfless guidance throughout the whole project. We sincerely thank the University of Edinburgh for its hospitality while VP was a visiting research student there for two months in 2016. Many thanks to Sverre Aarseth for his kind help with nbody6 and to Steve Shore for his perpetual support and thoughtful comments while finalising the paper. VP and this study was mainly supported by Charles University, grants GAUK186216 and SVV260441. LŠ also acknowledges support from the Czech Science Foundation through the project of Excellence No. 1437086G. We greatly appreciate access to the computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum, provided under the programme Projects of Large Research, Development, and Innovations Infrastructures (CESNET LM2015042). Large calculations were also performed on Tiger at the Astronomical Institute of Charles University. We thank the anonymous referee for useful comments.
References
 Aarseth, S. J. 1972, in Gravitational NBody Simulations, ed. M.Lecar (Cambridge, UK: Cambridge University Press), 88 [Google Scholar]
 Aarseth, S. J. 2003, Gravitational NBody Simulations (Cambridge, UK: Cambridge University Press) [CrossRef] [MathSciNet] [Google Scholar]
 Baumgardt, H., Heggie, D. C., Hut, P., & Makino, J. 2003, MNRAS, 341, 247 [NASA ADS] [CrossRef] [Google Scholar]
 Binney, J., & Tremaine, S. 1994, Galactic Dynamics (Princeton, NJ: Princeton University Press) [Google Scholar]
 Breen, P. G., & Heggie, D. C. 2012a, MNRAS, 425, 2493 [NASA ADS] [CrossRef] [Google Scholar]
 Breen, P. G., & Heggie, D. C. 2012b, MNRAS, 420, 309 [NASA ADS] [CrossRef] [Google Scholar]
 Casertano, S., & Hut, P. 1985, ApJ, 298, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1943, ApJ, 97, 255 [NASA ADS] [CrossRef] [Google Scholar]
 Cohn, H. 1980, ApJ, 242, 765 [NASA ADS] [CrossRef] [Google Scholar]
 Fujii, M. S., & Portegies Zwart, S. 2014, MNRAS, 439, 1003 [NASA ADS] [CrossRef] [Google Scholar]
 Giersz, M., & Heggie, D. C. 1994, MNRAS, 270, 298 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, J. 1987, ApJ, 313, 576 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, J., & Hut, P. 1993, ApJ, 403, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Heggie, D. C. 1975, MNRAS, 173, 729 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P. 1983, ApJ, 272, L29 [NASA ADS] [CrossRef] [Google Scholar]
 Hut, P. 1996, in Dynamical Evolution of Star Clusters: Confrontation of Theory and Observations, IAU Symposium, 174, 121 [NASA ADS] [Google Scholar]
 Larson, R. B. 1970, MNRAS, 150, 93 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D., & Eggleton, P. P. 1980, MNRAS, 191, 483 [NASA ADS] [Google Scholar]
 LyndenBell, D., & Wood, R. 1968, MNRAS, 138, 495 [NASA ADS] [Google Scholar]
 Makino, J. 1996a, in Dynamical Evolution of Star Clusters: Confrontation of Theory and Observations, IAU Symposium, 174, 151 [NASA ADS] [Google Scholar]
 Makino, J. 1996b, ApJ, 471, 796 [NASA ADS] [CrossRef] [Google Scholar]
 O’Leary, R. M., Stahler, S. W., & Ma, C.P. 2014, MNRAS, 444, 80 [NASA ADS] [CrossRef] [Google Scholar]
 Penston, M. V. 1969, MNRAS, 144, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Plummer, H. C. 1911, MNRAS, 71, 460 [NASA ADS] [CrossRef] [Google Scholar]
 Salpeter, E. E. 1955, ApJ, 121, 161 [NASA ADS] [CrossRef] [Google Scholar]
 Savitzky, A., & Golay, M. J. E. 1964, Anal. Chem., 36, 1627 [NASA ADS] [CrossRef] [Google Scholar]
 Spitzer Jr., L., & Hart, M. H. 1971, ApJ, 166, 483 [NASA ADS] [CrossRef] [Google Scholar]
 Storn, R., & Price, K. 1997, J. Global Optim., 11, 341 [CrossRef] [MathSciNet] [Google Scholar]
 Šubr, L. 2012, Astrophysics Source Code Library [record ascl:1206.007] [Google Scholar]
 Šubr, L., Kroupa, P., & Baumgardt, H. 2008, MNRAS, 385, 1673 [NASA ADS] [CrossRef] [Google Scholar]
 Takahashi, K. 1995, PASJ, 47, 561 [Google Scholar]
 Tanikawa, A., Hut, P., & Makino, J. 2012, New Astron., 17, 272 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Tables and figures
Correlation coefficient between the time of core collapse, t_{cc}, and the time of formation of the first binary with binding energy E_{bin} ≥ E_{lim} (i.e. t_{bin}, in our models; see also Figs. A.10 and A.12).
Fitted mean time of core collapse and the powerlaw index α of the core radius temporal evolution.
Mean logarithmic density gradients, a_{I − IV}, of our models at a given collapse.
Radii where the power law breaks, r_{I − III}, shown as logarithms to be consistent with the figures.
Fig. A.1. Lagrangian radii of model e10k. The mass fractions are on the right, the halfmass radius is plotted with a dashed line. Top: one arbitrary realisation. Smoothed curves are plotted in black, the original data are in grey. Bottom: average over all realisations. 

Open with DEXTER 
Fig. A.2. Same as in Fig. A.1 but for model e50k. 

Open with DEXTER 
Fig. A.3. Same as in Fig. A.1 but for model m20k. 

Open with DEXTER 
Fig. A.4. Same as in Fig. A.1 but for model m100k. 

Open with DEXTER 
Fig. A.5. Top: detail of the inner Lagrangian radii of one realisation of model e10k (left) and of model e50k (right). Circles correspond to the minima of smoothed radii, that is, set (8), and the dashed line is a powerlaw fit (11) through these data. Bottom: radial density profiles of this realisation at the given times in the range of radii plotted above. The dotted line demonstrates a fit by the triplebroken powerlaw function (15). The grey line shows an initial state of the system for comparison. 

Open with DEXTER 
Fig. A.6. Sequences of minima and binary star binding energy evolution in one realisation of model m20k(left) and of model m100k(right). Top: Lagrangian radii. Circles correspond to the minima of the radii, and the dashed lines are the powerlaw fits (11). Bottom: evolution of the binding energies of the dynamically formed binary stars. Each colour corresponds to one binary star family line, in which the individual components may be exchanged due to interactions with other stars. Black dashed vertical lines indicate the times of the homologous collapses, highlighted areas around them are the time intervals used for the evaluation of ΔE_{bin}. 

Open with DEXTER 
Fig. A.7. Same as in A.6 but for one realisation of model n100k. 

Open with DEXTER 
Fig. A.8. Distribution of the gain of binding energy, ΔE_{bin}, of the binary that hardened the most during the core collapse in model e10k. 

Open with DEXTER 
Fig. A.9. Distribution of the gain of binding energy, ΔE_{bin}, of the binary that hardened the most during a given homologous collapse in model m20k. 

Open with DEXTER 
Fig. A.10. Time of core collapse in model e10k versus the time of formation of the first binary of energy E_{bin} ≥ E_{lim} (indicated in each panel). The corresponding correlation coefficients are listed in Table A.1. 

Open with DEXTER 
Fig. A.11. Same as in Fig. A.10 but for model e50k. 

Open with DEXTER 
Fig. A.12. Same as in Fig. A.10 but for model m20k. The corresponding correlation coefficients are listed in Table A.1. 

Open with DEXTER 
Fig. A.13. Same as in Fig. A.10 but for model m100k. 

Open with DEXTER 
All Tables
Correlation coefficient between the time of core collapse, t_{cc}, and the time of formation of the first binary with binding energy E_{bin} ≥ E_{lim} (i.e. t_{bin}, in our models; see also Figs. A.10 and A.12).
Fitted mean time of core collapse and the powerlaw index α of the core radius temporal evolution.
Mean logarithmic density gradients, a_{I − IV}, of our models at a given collapse.
Radii where the power law breaks, r_{I − III}, shown as logarithms to be consistent with the figures.
All Figures
Fig. 1. Schematic plot of the radial density profile in an Nbody cluster during the core collapse with four different values of the logarithmic density gradient, a_{I − IV}. The dotted line has a slope equal to α and shows the asymptotic solution of LyndenBell & Eggleton (1980). The break radius r_{I} is identified with the core radius, r_{III} roughly corresponds to the halfmass radius, and r_{k} is the cluster radius. The notation is the same as in Eq. (15) that we used for the fitting of our numerical models. 

Open with DEXTER  
In the text 
Fig. A.1. Lagrangian radii of model e10k. The mass fractions are on the right, the halfmass radius is plotted with a dashed line. Top: one arbitrary realisation. Smoothed curves are plotted in black, the original data are in grey. Bottom: average over all realisations. 

Open with DEXTER  
In the text 
Fig. A.2. Same as in Fig. A.1 but for model e50k. 

Open with DEXTER  
In the text 
Fig. A.3. Same as in Fig. A.1 but for model m20k. 

Open with DEXTER  
In the text 
Fig. A.4. Same as in Fig. A.1 but for model m100k. 

Open with DEXTER  
In the text 
Fig. A.5. Top: detail of the inner Lagrangian radii of one realisation of model e10k (left) and of model e50k (right). Circles correspond to the minima of smoothed radii, that is, set (8), and the dashed line is a powerlaw fit (11) through these data. Bottom: radial density profiles of this realisation at the given times in the range of radii plotted above. The dotted line demonstrates a fit by the triplebroken powerlaw function (15). The grey line shows an initial state of the system for comparison. 

Open with DEXTER  
In the text 
Fig. A.6. Sequences of minima and binary star binding energy evolution in one realisation of model m20k(left) and of model m100k(right). Top: Lagrangian radii. Circles correspond to the minima of the radii, and the dashed lines are the powerlaw fits (11). Bottom: evolution of the binding energies of the dynamically formed binary stars. Each colour corresponds to one binary star family line, in which the individual components may be exchanged due to interactions with other stars. Black dashed vertical lines indicate the times of the homologous collapses, highlighted areas around them are the time intervals used for the evaluation of ΔE_{bin}. 

Open with DEXTER  
In the text 
Fig. A.7. Same as in A.6 but for one realisation of model n100k. 

Open with DEXTER  
In the text 
Fig. A.8. Distribution of the gain of binding energy, ΔE_{bin}, of the binary that hardened the most during the core collapse in model e10k. 

Open with DEXTER  
In the text 
Fig. A.9. Distribution of the gain of binding energy, ΔE_{bin}, of the binary that hardened the most during a given homologous collapse in model m20k. 

Open with DEXTER  
In the text 
Fig. A.10. Time of core collapse in model e10k versus the time of formation of the first binary of energy E_{bin} ≥ E_{lim} (indicated in each panel). The corresponding correlation coefficients are listed in Table A.1. 

Open with DEXTER  
In the text 
Fig. A.11. Same as in Fig. A.10 but for model e50k. 

Open with DEXTER  
In the text 
Fig. A.12. Same as in Fig. A.10 but for model m20k. The corresponding correlation coefficients are listed in Table A.1. 

Open with DEXTER  
In the text 
Fig. A.13. Same as in Fig. A.10 but for model m100k. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext 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 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.