CATCH ME IF YOU CAN: IS THERE A "RUNAWAY-MASS" BLACK HOLE IN THE ORION NEBULA CLUSTER?

The canonical (best-fit) model of the ONC has an initial mass of 5400 M_☉, half of which is in the form of stars, while gas accounts for the other half. We considered a Kroupa (2001) IMF, which, for the given cluster star mass, predicts ≈50 OB stars to have been formed in the ONC (Pflamm-Altenburg & Kroupa 2006). In order to avoid additional statistical noise in the evolution of the number of OB stars, we used identical samples of stellar masses for all numerical realizations of the model. The initial half-mass radius of the cluster is r_h ≈ 0.11 pc, which is compatible with the initial or birth radius–mass relation inferred for star clusters ranging in mass up to globular clusters (Marks & Kroupa 2012). Positions and velocities are generated in a mass segregated state according to the prescription of Šubr et al. (2008) with mass segregation index S = 0.4. The algorithm by definition places massive stars in the cluster core and, therefore, the initial half-mass radius of the OB stars is ≈0.05 pc. In order not to place the light gas particles at the cluster outskirts, we generated initial positions for stellar and gas particles separately, i.e., we had a constant gas to star mass ratio throughout the whole cluster at T = 0. All OB stars were set to be members of primordial binaries with a secondary mass M_s ⩾ 1 M_☉ with the pairing algorithm being biased toward assigning a massive secondary to a massive primary. More specifically, the algorithm first sorts the stars from the most massive to the lightest one. The most massive star from the set is taken as the primary. The secondary star index, id, in the ordered set is generated as a random number with the probability density ∝id^−β and max (id) corresponding to a certain mass limit, M_{s, min}. The two stars are removed from the set and the whole procedure is repeated until stars with M_⋆ ⩾ M_{p, min} remain. Typically, we used as the minimal mass of the primary M_{p, min} = 5 M_☉, as the minimal mass of the secondary M_{s, min} = 1 M_☉ and the pairing algorithm index β = 40. For the assigned binaries we used a semimajor axis distribution according to Öpik's law n(a)∝a⁻¹ (Kobulnicky & Fryer 2007) and a thermal distribution of eccentricities, n(e)∝e. Except for the possibility of being a secondary member of a binary containing an OB star, we did not generate binaries of low-mass stars. This limitation comes from the requirement of numerical effectivity and it is not likely to affect our results substantially.

Let us note that the canonical model presented above fits into the range of the expanding class of models of the ONC discussed by Kroupa (2000). Other possible initial conditions would be fractal models that collapse (e.g., Allison et al. 2010). These require star formation to be synchronized across the pre-cluster cloud core to much shorter than the dynamical time though.

Most of the features of the evolution of a star cluster with strong few-body relaxation and gas expulsion can be demonstrated with the canonical model. In order to distinguish systematic effects from rather large fluctuations of cluster parameters, we present quantities averaged over 100 realizations of the model. The upper panel of Figure 1 shows the temporal evolution of the half-mass radius of the stellar component, r_h. During the initial phase, the cluster expands due to two-body relaxation. Gas expulsion that starts at T_ex = 0.5 Myr removes all gas from the cluster within a few hundred thousand years. From the point of view of the stars, this event leads to an abrupt decrease of their potential energy. Consequently, they can reach larger distances from the cluster center, which manifests itself as an accelerated growth of r_h at 0.5 Myr ≲ T ≲ 1 Myr. Initially weakly bound stars became unbound due to the gas expulsion and they escape from the cluster. Hence, reduction of the cluster mass also accelerates—see the lower panel of Figure 1 where we plot the mass of the star cluster. We define M_c as the sum of the mass of the stars within a sphere of radius r_lim from the cluster center. We set r_lim = 3 pc in order to match typical observations of the ONC, which usually count stars within a projected distance of ≈3 pc from the Trapezium. Note that the decrease of M_c is accelerated with ≈0.5 Myr delay after the time of gas expulsion. This is due to the unbound stars taking some time to reach r_lim. Besides the escape of stars due to the gas expulsion, there is also continuous stellar mass loss due to the few-body relaxation, which is capable of accelerating stars above the escape velocity. From T ≈ 1 Myr onward the cluster expansion is again dominated by relaxational processes driving a revirialization (Kroupa et al. 2001).

Zoom In Zoom Out Reset image size

Our model indicates that the initial mass of the ONC must have been at least 20% larger than it is now, after a few million years of dynamical evolution. We also see that, mainly due to the effect of gas expulsion, the ONC must have been more compact by a factor ≳ 5 at the time of its birth than it is now. Hence, the initial relaxation time was much shorter than what we would infer from present-day observations. Consequently, close few-body interactions must have been more frequent in the past. In particular, we suggest that scattering of single stars on binaries is a process that has played an important role in the evolution of the ONC. This kind of interaction has been studied intensively in the literature (e.g., Heggie & Hut 2003 and references therein). A rich variety of outcomes can be obtained, depending on the initial conditions. Nevertheless, some general results can be formulated. In particular, if the binary is hard, i.e., its orbital velocity is larger than the impact velocity of the third star, it typically shrinks, losing its energy, which is transferred to the accelerated single star. Therefore, we expect two processes to happen in correlation: (1) ejections of high-velocity stars and (2) physical collisions of massive binaries, which lead to the formation of more massive objects.

The probability of collisions increases with the mass of the interacting stars. Therefore, numerical models by different authors often show a "runaway" process when most of the collisions involve the most massive object which grows continuously (e.g., Portegies-Zwart et al. 2004). The same is true also for our model of the ONC. The upper panel of Figure 2 shows the growth of the mass of the most massive object, M_•, due to stellar collisions; the number of OB stars left in the cluster is plotted in the lower panel. Approximately, one third of the missing OB stars have disappeared due to merging, while the other two thirds have been ejected from the cluster with velocities >10 km s⁻¹. The merging tree varies significantly among individual realizations of the particular model. In general, more than one merging star is formed during the initial ≈0.5 Myr. Among other processes, this is due to the merging of several primordial OB binaries. At later stages, usually one runaway-mass object dominates the merging process.

Zoom In Zoom Out Reset image size

In the Appendix we provide an approximate description of the rate of stellar ejections due to the scattering of single stars on massive binaries. It gives an estimate of the decay time of the number of the OB stars of τ ≳ 5 Myr for the canonical model. If we assume that the decay rate is linearly proportional to the number of OB stars, we expect roughly an exponential decay of N_OB, i.e., it should reach half of its initial value at ∼3.5 Myr, which is in good agreement with the numerical results.

A mean number of OB stars of ≈28 remain at T = 2.5 Myr. This considerably exceeds the number of OB systems observed in the ONC—the study of Hillenbrand (1997) reports only 10 stars heavier than 5 M_☉. However, the real number of massive stars may be somewhat larger due to their "hiding" in OB binary systems. Our models always end up with several OB stars having an OB companion and, therefore, a better agreement with the observations is achieved if we count the OB systems as indicated with the dashed line in Figure 2. Moreover, due to the stochastic nature of dynamical cluster evolution, different realizations of the model end up with quite different values of N_OB and M_•. The final state (T = 2.5 Myr) of all realizations of the canonical model is shown in Figure 3. Although the distribution of points in the N_OB–M_• space is rather noisy, we can deduce an anticorrelation between these two quantities (their linear correlation coefficient is −0.42). In particular, all realizations that end up with N_OB ⩽ 20 lead to the formation of a runaway-mass star of mass M_• > 100 M_☉. In other words, the underabundance of the high-mass stars in the ONC not only indicates a period of prominent two-body relaxation in the past, but also the merging formation of a massive object that may represent an important footprint of the cluster's history.

Zoom In Zoom Out Reset image size

Interestingly, the term "runaway" has a secondary meaning in some cases. Despite of its high mass, in several realizations the merging object has been ejected out of the cluster with velocity exceeding 10 km s⁻¹. These cases can be identified due to differences of the upper and lower panels of Figure 3. The escape of the most massive body beyond the 3 pc boundary is also responsible for the drops of the dashed line in the plot of the temporal evolution of M_• in Figure 2, which shows the mean mass of the most massive object within 3 pc from the cluster center.

As mentioned above, the canonical model indicates an anticorrelation between the final number of the OB stars and the final mass of the runaway-mass star. In Figure 4 we plot the mean mass of the most massive object versus the mean number of remaining OB stars for all models listed in Table 1. Now, each point represents an average over several tens of realizations of the particular model and an N_OB–M_• anticorrelation becomes evident. We attribute this relation to the fact that both mechanisms of OB star removal, i.e., physical collisions and ejections, are driven by a common underlying process of close three-body interactions (see the Appendix). Naturally, their importance grows with increasing stellar density. Hence, the models that stay more compact during the course of their evolution are located in the top-left corner of the graph. They better fit the observations from the point of view of the number of OB stars, however, either their half-mass radius or total stellar mass at T = 2.5 Myr is too small, i.e., not consistent with the observations.

Zoom In Zoom Out Reset image size

In most of our models, the gas particles were blown out of the cluster due to a repulsive external force $\propto {\bm r}/r^3$ centered on the cluster core. The strength of the force was set such that the gas particles were accelerated to velocities ≈8 km s⁻¹ on the timescale of ≈0.1 Myr. This approach is the most realistic one, keeping the continuity equation fulfilled, but, at the same time, it is numerically the most expensive. The alternative method with instantaneous removal of the gas particles at T_ex leads to similar results. The clusters integrated with this method typically result in somewhat lower total mass and larger half-mass radius at the final time. Let us note, for completeness, that in the literature, another method that mimics gas removal in numerical models of star clusters is also used (e.g., Geyer & Burkert 2001; Kroupa et al. 2001; Baumgardt & Klessen 2011). It is based on the time-variable smooth external potential that represents the weakening gravitational potential of the gas. It has been found that the two approaches (i.e., modeling gas via particles versus external potential) lead to nearly indistinguishable results (Geyer & Burkert 2001). In order to check for a possible bias of our gas particles approach that may stem from two-body relaxational processes, we have integrated two additional models with very low masses for the gas particles. Due to the large number of particles in these integrations, we have followed their evolution only up to the time of the gas expulsion, T_ex = 0.5 Myr. Outcomes of these models in terms of the mean number of remaining OB stars, N_OB, and the mass of the runaway-mass star, M_•, are presented in Table 2. Besides just giving the mean values of N_OB and M_•, we have also performed a Kolmogorov–Smirnov test comparing the sets of N_OB and M_• coming from these models with the canonical one. As can be seen, there is no apparent dependence of the results on m_gas ranging from 0.03 to 0.4 M_☉. Hence, we conclude that the two-body relaxation due to the gas particles does not play a significant role in our simulations.

Table 2. Comparison of Models with Different Masses of Gas Particles

mg	NOB	M•	Nrun	$\!\! P_\mathrm{KS,N_{OB}} \!\!$	$\!\! P_\mathrm{KS,M_\bullet } \!\!\!$	Comment
0.40	39.7	125.1	40	95%	42%	Model 4
0.20	39.8	124.9	100	⋅⋅⋅	⋅⋅⋅	Canonical model
0.05	38.2	128.5	14	82%	87%
0.03	39.9	126.7	14	32%	75%

Notes. Mean values of the number of OB stars within the radius of 3 pc from the cluster center, N_OB, and mass of the runaway-mass star, M_•, for models with different values of m_g. All other parameters of the models, including the total mass of gas, are identical to the canonical model. P_{KS, *} is the Kolmogorov–Smirnov probability that the null hypothesis (assuming the particular data set comes from the same distribution as the canonical one) is valid. Note that a value of P_{KS, *} ≈ 5% is usually considered to be a limit below which the null hypothesis should be rejected.

Download table as:  ASCIITypeset image

The evolution of the radiation pressure in the real ONC has been definitely more complicated. According to observations (Huff & Stahler 2006), the star formation in the ONC was continuous with the most massive stars starting to be formed no more than 2 Myr ago. Setting T_ex in our models to either 0.5 or 0.7 Myr appears to have only a marginal impact on the final state of the cluster. Hence, we assume that a more complicated temporal prescription of the gas removal would not affect our results considerably.

Under the assumption of an invariant canonical IMF and a small initial cluster radius, the formation of the runaway-mass object via stellar collisions appears to be an inevitable process that, together with the high-velocity ejections, decreases the number of OB stars in a dense star cluster. Its possible detection would definitely strongly support the scenario of the ONC history presented in this paper.

The current state of the runaway-mass star is, however, not clear. The massive star would definitely have undergone a fast internal evolution. While successive merging events may lead to very fast mass growth, stellar winds act in the opposite way. The state of the runaway-mass object after ≈2 Myr of evolution depends on the (dis)balance of these two processes. In the literature, quite different predictions on this subject can be found. Glebbeek et al. (2009) suggest that winds should work in a self-regulatory manner, limiting the maximum mass and, consequently the star's lifetime. Other authors (e.g., Suzuki et al. 2007; Pauldrach et al. 2011) state that stellar winds of the runaway-mass star will not be able to terminate its growth and it may then collapse to a massive black hole. There is no observational evidence for a star with a mass above 50 M_☉ in the ONC that could be interpreted as a runaway-mass object. The most massive stellar member of the ONC, Θ¹C, is a binary with a total mass of about 45 M_☉ and mass ratio ≈0.25 (Kraus et al. 2009), i.e., the primary has mass ≲ 35 M_☉. It does not exhibit a stellar wind strong enough to reduce its mass by several tens of M_☉. Hence, our scenario assumes a low efficiency of stellar winds, i.e., continuous growth of the runaway-mass star, which forms a massive black hole at the end of its lifetime. There is no general consensus regarding whether the collapse of a star more massive than 150 M_☉ will be followed by an ejection of most of its mass. As there is no evidence for a supernova remnant in the ONC, we assume most of the mass of the runaway-mass star, if it is present, to have collapsed directly into a black hole.

The putative black hole in the ONC would be detectable only through an interaction with its environment. One possibility is an accretion of surrounding gas, which could be detectable through X-ray radiation. The main source of the gas in the Trapezium region is stellar winds from the remaining OB stars. We estimate that they produce at most 10⁻⁵ M_☉ yr⁻¹ within the central 0.2 pc. However, the black hole is likely to capture only the gas that falls within its Bondi radius ≈100 AU, i.e., the accretion rate could be ≲ 10⁻¹⁵ M_☉yr⁻¹. Such a highly sub-Eddington accretion will not be detectable in the region confused with several X-ray young stellar sources.

A better chance for detecting the black hole would be given if it is a member of a binary. This indeed appears to be the case for about two thirds of the realizations of the canonical model. However, the typical separation of the black hole and the secondary appears to be from several tens to hundreds of AU (see Figure 5). Even with the relatively large eccentricities achieved, in no case does the secondary star fill its Roche radius at the pericentre of its orbit, i.e., we do not expect the black hole to be a member of a mass transferring system. Still there would be a chance for a detectable accretion provided the secondary star is massive (M_⋆ ≳ 10 M_☉). In such a case, its stellar wind may produce short periods of enhanced activity during pericentre passages, but the corresponding luminosity is rather uncertain. Let us consider the well known Cyg-X1 system as a kind of template. Having a black hole mass of ≈10 M_☉ and separation from the secondary star of ≈0.2 AU, it has a luminosity ≈2 × 10³⁷ erg s⁻¹ (Wilms et al. 2006), which corresponds to ≈10⁻⁴ of the rest-mass energy of the winds emitted by the secondary star (≈10⁻⁶ M_☉ yr⁻¹). Assuming the same effectivity in conversion of the stellar winds into radiation, we estimate the luminosity to be similar to the Cyg-X1 system, i.e., ≈10⁴ L_☉, for separations of the order of 1 AU. Considering the density of the stellar wind to decrease with the square of the distance, we obtain a rough estimate of the peak luminosity L ≈ 10⁴(r_p/1 AU)⁻² L_☉. According to the mass of the putative black hole, the maximum of the emitted radiation is expected to lie in the X-ray band. The accretion events should be recurrent with a period ranging from years to hundreds of years. A shorter period generally implies a smaller r_p, i.e., larger values of the peak luminosity.

Zoom In Zoom Out Reset image size

Another way of an indirect detection of the massive black hole in the ONC lies in velocity dispersion measurements. It has been stated by several authors that the core of the ONC (r ≲ 0.25 pc) is dynamically hot with a velocity dispersion ≳ 4 km s⁻¹ (e.g., Jones & Walker 1988; van Altena et al. 1988; Tobin et al. 2009). The missing stellar mass required for virial equilibrium in the innermost region of the ONC was estimated to be ≳ 2000 M_☉ which exceeds the observed mass ≈200 M_☉ by an order of magnitude. Kroupa et al. (2001) demonstrate that the globally super-virial velocity dispersion is readily obtained if the ONC is expanding now after expulsion of its residual gas, but some of the missing mass in the inner region could be attributed to the invisible remnant of the runaway-mass star.

As we show in Figure 6, our canonical model gives a velocity dispersion ≳ 4 km s⁻¹ in the inner 0.25 pc, which is in accord with the observational data. In order to inspect the influence of the mass of the runaway-mass object on the velocity dispersion, we have divided the 100 realizations of the canonical model into two groups of 50 according to the mass of the most massive object that remained in the cluster. It appears that the clusters with M_• ≳ 120 M_☉ have a velocity dispersion within 0.25 pc somewhat larger than those with M_• ≲ 120 M_☉, nevertheless, even within the latter group we have σ ≳ 3.5 km s⁻¹, i.e., an apparently super-virial core region. Together with the fact that none of our numerical experiments led to M_• ≳ 2000 M_☉, virtually required by the condition of virial equilibrium, this indicates that the large velocity dispersion can be only partly due to the hidden mass of the runaway-mass object. Detailed analysis of our models shows that there are two other reasons for having large velocity dispersion in the core. First, according to our models, the ONC is a post-core-collapse cluster with centrally peaked density. This leads to a deep gravitational potential well that allows the central velocity dispersion to be larger than what would be expected, e.g., for the Plummer profile. Second, in our models and perhaps also in the real ONC, there are likely to be present several unidentified wide binaries with orbital velocities exceeding 10 km s⁻¹, which affect the velocity dispersion.

Zoom In Zoom Out Reset image size

Figure 6 indicates that the presence of the dark massive body influences remarkably (i.e., by more than 1 km s⁻¹) the velocity dispersion in a small region of a radius of ≈0.1 pc. This is in accord with an analytical estimate of the influence radius of the black hole, r_h = GM_•/σ², which for M_• = 100 M_☉ and σ = 3 km s⁻¹ gives r_h ≈ 0.05 pc. Due to the small number count of stars in such a small region, standard statistical methods do not represent a robust tool for determining the presence of a dark massive body. Nevertheless, having this caution in mind, let us discuss the velocity dispersion of the compact core of the ONC, the four Trapezium stars. According to the published observational data (see Table 3), the line-of-sight velocity dispersion of the Trapezium is σ_los ≈ 4.6 km s⁻¹, which implies a three-dimensional velocity dispersion σ ≈ 7.9 km s⁻¹. Considering a sphere of radius r_T ≈ 0.025 pc covering the four Trapezium stars, we obtain an estimate of the enclosed mass required for the system to be in virial equilibrium: M_bind ≈ r_T σ²/G ≈ 350 M_☉. Taking into account the stellar mass of the Trapezium (≈90 M_☉), we find it to be an apparently dynamically very hot system. Our models suggest that the virial equilibrium state can be achieved by the presence of the runaway-mass object.

Yet another, and statistically better founded, argument for a hidden massive body comes from the comparison of the kinematical state of the Trapezium and similar subsystems found in our numerical models. We have implemented an algorithm similar to that used by Allison & Goodwin (2011) for detection of Trapezium-like systems. In particular, for each OB star we determined its three nearest OB neighbors (either single OB stars, or bound systems containing at least one OB star). The set of stars was considered Trapezium-like if all members had the same OB neighbors. In order to obtain a statistically significant sample, we examined the last 20 snapshots (covering the period from ≈1.5 Myr to 2.5 Myr) of 100 independent realizations of the canonical model. We found Trapezium-like systems in a third of cluster snapshots. The mean value of the velocity dispersion of these systems was found to be ≈7.5 km s⁻¹. The runaway-mass star was excluded from the search algorithm. Furthermore, we have distinguished systems according to their center of mass distance to the runaway-mass object. In approximately two thirds of the Trapezium-like systems, the runaway-mass object was found within a distance smaller than max (Δr_ij), with Δr_ij denoting separations of individual members. The mean velocity dispersion of the Trapezium-like system of this subset was ≈9.0 km s⁻¹, while it was only ≈4.3 km s⁻¹ for the remaining cases.

Finally, let us remark that kinematical evidence of a massive black hole in the ONC may also lie in the ≈70% probability that it is a member of a binary (see Section 3.3). The typical velocity of the secondary star should be >10 km s⁻¹, i.e., considerably exceeding the observed central velocity dispersion.

The frequency of scattering events of a star with velocity v_c on binaries with number density n_B and semimajor axis a is

$$\begin{equation} \nu = \Sigma \,n_\mathrm{B}\,v_\mathrm{c} = \frac{2\pi G M_\mathrm{B} n_\mathrm{B} a}{v_\mathrm{c}}, \end{equation} \tag{ A7 }$$

where we assumed the binary cross-section to be determined by the gravitational focusing, i.e., Σ ≈ 2πGM_Ba/v²_c. The mean frequency of ejections can be obtained via integration of Equation (A7) weighted by the distribution function (A5):

$$\begin{equation} \bar{\nu }_\mathrm{e} = \int _{M_\mathrm{min}}^{M_\mathrm{max}} {d}M_\mathrm{B}\int _{a_\mathrm{min}}^{a_\mathrm{lim}} {d}a\, \nu (a,M_\mathrm{B})\, n(a, e, M_\mathrm{B}). \end{equation} \tag{ A8 }$$

The upper limit of the semimajor axis, a_lim, has to be set such that only interactions that lead to acceleration above the escape velocity from the cluster are considered, i.e.,

$$\begin{equation} v_\mathrm{acc}(a_\mathrm{lim}) = v_\mathrm{esc} \approx \sqrt{\frac{2G M_\mathrm{c}}{r_\mathrm{c}}}\;, \end{equation} \tag{ A9 }$$

where r_c is a characteristic radius of the cluster. Combining Equations (A9) and (A4) gives a_lim ≈ r_cM_B/M_c. For the canonical model with initial M_c = 5400 M_☉ and r_c = 0.11 pc and for 5 M_☉ ⩽ M_B ⩽ 100 M_☉ we obtain a_lim ≳ 100 AU, which is the upper limit of the semimajor axis distribution used in our numerical integrations. Hence, let us for simplicity assume a_lim = a_max = const., which yields

$$\begin{equation} \bar{\nu }_\mathrm{e} = \left\lbrace \begin{array}{@{}lll@{}} \displaystyle \frac{2\pi G n_\mathrm{B}}{v_\mathrm{c}}\, \frac{a_\mathrm{max} - a_\mathrm{min}}{\ln (a_\mathrm{max}/a_\mathrm{min})}\, \frac{M_\mathrm{max} - M_\mathrm{min}}{\ln (M_\mathrm{max}/M_\mathrm{min})} & {\rm for} & \alpha =1, \\\ms \displaystyle\frac{2\pi G n_\mathrm{B}}{v_\mathrm{c}}\, \frac{a_\mathrm{max} - a_\mathrm{min}}{\ln (a_\mathrm{max}/a_\mathrm{min})}\, \frac{1 - \alpha }{2 - \alpha }\, \frac{M_\mathrm{max}^{2-\alpha } - M_\mathrm{min}^{2-\alpha }}{M_\mathrm{max}^{1-\alpha } - M_\mathrm{min}^{1-\alpha }} & {\rm for} & \alpha \ne 1. \end{array} \right. \end{equation} \tag{ A10 }$$

Considering an initial OB binary density n_B ≈ 10⁵ pc⁻³ and velocity dispersion v_c ≈ 10 km s⁻¹, a ∈ 〈0.1 AU, 100 AU〉, and M_B ∈ 〈5 M_☉, 100 M_☉〉 with a Salpeter mass function (α = 2.35) we obtain $\bar{\nu }_\mathrm{e} \approx 0.1 \,\mathrm{Myr}^{-1}$. Hence, the mean time for an OB star to undergo scattering on a massive binary that leads to its ejection from the cluster is $\tau _\mathrm{e} \equiv 1 / \bar{\nu }_\mathrm{e} \approx 10\,\mathrm{Myr}$ for the canonical model.

Unlike in the case of stellar ejections, we assume that stellar collisions are caused by binary shrinking due to successive three-body interactions. We further assume all the massive binaries to be hard, i.e., their interaction with the third body leads to the acceleration of the impact star and growth of the binding energy of the binary. The frequency of the events of scattering of a star of mass M_⋆ on the massive binary is ν' = Σv_cn(M_⋆)dM_⋆, where n(M_⋆)dM_⋆ is the number of stars of given mass per unit volume. The rate of the energy transfer from the binary to the stars of mass in 〈M_⋆, M_⋆ + dM_⋆〉 is

$$\begin{equation} \frac{{d}E}{{d}t} \approx \Delta E_\star \nu ^\prime = \frac{2\pi G^2 M_\mathrm{B}^2 M_\star }{v_\mathrm{c}}\, n(M_\star){d}M_\star \;. \end{equation} \tag{ A11 }$$

The net rate of energy transfer to stars within the whole mass spectrum gives

$$\begin{equation} \frac{{d}E}{{d}t} \approx \int \frac{2\pi G^2 M_\mathrm{B}^2 M_\star }{v_\mathrm{c}} n(M_\star) {d}M_\star = \frac{2\pi G^2 M_\mathrm{B}^2 \rho _\star }{v_\mathrm{c}} \;. \end{equation} \tag{ A12 }$$

The stars collide when the binary separation at pericentre becomes less than the sum of the stellar radii. In terms of binding energy and with approximation R_⋆(M₁) + R_⋆(M₂) ≈ R_⋆(M_B), the condition for collision reads

$$\begin{equation} E_\mathrm{coll} \approx \frac{1-e}{8} \frac{G M_\mathrm{B}^2}{R_\star (M_\mathrm{B})}\;. \end{equation} \tag{ A13 }$$

The time required to grow the binding energy from the initial value E₀ ≈ (1/8)GM²_B/a₀ to E_coll is

$$\begin{equation} t_\mathrm{coll} \approx (E_0 - E_\mathrm{coll}) \left| \frac{{d}E}{{d}t} \right|^{-1} \approx \frac{v_\mathrm{c}}{8\pi G \rho _\star } \left(\frac{1-e}{R_\star (M_\mathrm{B})} - \frac{1}{a_0} \right)\;. \end{equation} \tag{ A14 }$$

Integration of t_coll over the whole parameter space of binaries with the distribution function (A5) gives the characteristic time of stellar collisions:

$$\begin{eqnarray} \tau _\mathrm{coll} & \approx & \frac{v_\mathrm{c}}{16\pi G \rho _\star } \int _{a_\mathrm{min}}^{a_\mathrm{max}} \! {d}a \int _{M_\mathrm{min}}^{M_\mathrm{max}} \! {d}M_\mathrm{B} \int _0^{1-R_\star (M_\mathrm{B})/a} \! {d}e \left(\frac{1-e}{R_\star (M_\mathrm{B})} - \frac{1}{a_0} \right) n(a, e, M_\mathrm{B}) \nonumber \\ & \approx & \frac{v_\mathrm{c}}{16\pi G \rho _\star } \int _{a_\mathrm{min}}^{a_\mathrm{max}} \! {d}a \int _{M_\mathrm{min}}^{M_\mathrm{max}} \! {d}M_\mathrm{B} \, A\,a^{-1} M_\mathrm{B}^{-\alpha } \left(\frac{1}{3R_\mathrm{B}} - \frac{1}{a} + \frac{R_\mathrm{B}}{a^2} - \frac{R_\mathrm{B}^3}{3a^3} \right)\;, \end{eqnarray} \tag{ A15 }$$

where we denoted R_B ≡ R_⋆(M_B); the upper limit on the eccentricity is set such that the binaries are not collisional initially. Assuming R_B ≫ a, we omit the two least significant terms in Equation (A15) and perform the integration over a:

$$\begin{equation} \tau _\mathrm{coll} \approx \frac{A v_\mathrm{c}}{16\pi G \rho _\star } \int _{M_\mathrm{min}}^{M_\mathrm{max}} \! {d}M_\mathrm{B}\, M_\mathrm{B}^{-\alpha } \left[ \frac{1}{3R_\mathrm{B}}\ln \left(\frac{a_\mathrm{max}}{a_\mathrm{min}} \right) + \frac{1}{a_\mathrm{max}} - \frac{1}{a_\mathrm{min}} \right]\;. \end{equation} \tag{ A16 }$$

Finally, using the notation $\tilde{M} \equiv M/M_\odot$ we obtain for $R_\star (M) = R_\odot \,\tilde{M}^{0.8}$:

$$\begin{equation} \tau _\mathrm{coll} \approx \left\lbrace \begin{array}{@{}lll@{}} \displaystyle \frac{v_\mathrm{c}}{16\pi G \rho _\star R_\odot } \left[ \frac{5}{12}\left(\tilde{M}_\mathrm{min}^{-0.8} - \tilde{M}_\mathrm{max}^{-0.8} \right) \ln ^{-1} \! \left(\frac{M_\mathrm{max}}{M_\mathrm{min}} \right) + \left(\frac{R_\odot }{a_\mathrm{max}} - \frac{R_\odot }{a_\mathrm{min}} \right) \ln ^{-1} \! \left(\frac{a_\mathrm{max}}{a_\mathrm{min}} \right) \right] & {\rm for} & \alpha = 1, \\\ms \displaystyle \frac{v_\mathrm{c}}{16\pi G \rho _\star R_\odot } \left[ \frac{1-\alpha }{0.6-3\alpha } \frac{\tilde{M}_\mathrm{max}^{0.2-\alpha } - \tilde{M}_\mathrm{min}^{0.2-\alpha }}{\tilde{M}_\mathrm{max}^{1-\alpha } - \tilde{M}_\mathrm{min}^{1-\alpha }} + \left(\frac{R_\odot }{a_\mathrm{max}} - \frac{R_\odot }{a_\mathrm{min}} \right) \ln ^{-1} \! \left(\frac{a_\mathrm{max}}{a_\mathrm{min}} \right) \right] & {\rm for} & \alpha \ne 1. \end{array} \right. \end{equation} \tag{ A17 }$$

The canonical model has an initial central density ρ_⋆ ≈ 2 × 10⁶ M_☉ pc⁻³, which gives τ_coll ≈ 40 Myr, i.e., we estimate stellar collisions to be roughly a factor of four less efficient process for OB star removal than three-body scattering. This is somewhat less than what our numerical experiment shows. This discrepancy can be due to several reasons. Most important is probably the fact that stellar collisions may occur sooner due to perturbations that lead to a strong growth of orbital eccentricity. Indeed, the numerical experiments show that most of the collisions occur with e ≳ 0.99. Hence, the above given derivation can only serve as an order of magnitude estimate.

Altogether, the three-body interactions are expected to decrease the number of OB stars in the canonical model of the ONC by a factor of two on a timescale of ≈5 Myr, which is in accord with the results of the numerical experiment.