DEPENDENCY OF DYNAMICAL EJECTIONS OF O STARS ON THE MASSES OF VERY YOUNG STAR CLUSTERS

The majority of O stars are found in young star clusters or OB associations. It has been suggested that a small fraction of massive stars may form in isolation as a result of stochastic sampling of the stellar and cluster mass function (see de Wit et al. 2005; Parker & Goodwin 2007; Oey & Lamb 2012). However, in the case of the Galactic field O stars within 2–3 kpc from the Sun, essentially all of them can be traced back to young star clusters or associations (Schilbach & Röser 2008), suggesting the field O stars also formed in clusters (Gvaramadze et al. 2012). Some observational studies on O stars in the Large and Small Magellanic Clouds have claimed there to be evidence for isolated massive star formation (Bressert et al. 2012; Oey et al. 2013), such as e.g., apparent isolation, absence of a bow-shock, and circular shape of an H ii region. These claims cannot, however, be conclusive since they can be explained with slowly moving runaways or former members of clusters which have already dissolved (Weidner et al. 2007). For example, only 30–40% of the Galactic OB runaways are found to have a bow-shock (van Buren et al. 1995; Huthoff & Kaper 2002), and, furthermore, bow-shock formation depends on the physical condition of the ambient medium through which a runaway travels (e.g., temperature and density of the ambient medium, Huthoff & Kaper 2002). Mackey et al. (2013) showed that runaways can have a circular H ii region in the projected Hα emission. Furthermore, in the case of the 30 Dor region Bressert et al. considered, there are many young star clusters in the region. If the 30 Dor region is modelled by a population of clusters, how many O stars would be ejected or removed from their clusters in 10 Myr? This has not been investigated. Any conclusions on the issue of isolated O star formation are, for the time being, thus premature. Finally, measuring stellar proper motions in the Large or Small Magellanic Clouds is unfeasible such that claims for O stars formed in isolation become next to impossible to confirm or reject. With Gaia - , however, it may become possible to put constraints on this if proper motions of massive stars can be measured to a precision of a few km s⁻¹ in the Magellanic Clouds. Given the above results on the well studied ensemble of O stars within 2–3 kpc around the Sun it appears physically more plausible that all O stars form in embedded clusters.

There are two mechanisms to expel O stars from a star cluster. One is the binary supernova scenario (Blaauw 1961; Portegies Zwart 2000; Eldridge et al. 2011) in which a star in a binary obtains a high velocity after the supernova explosion of the initially more massive companion star. The other is dynamical ejection through strong close encounters with other massive stars in the cluster core (Poveda et al. 1967; Leonard & Duncan 1990; Fujii & Portegies Zwart 2011; Banerjee et al. 2012; Perets & Šubr 2012). Of the candidate isolated O stars cataloged by Schilbach & Röser (2008), 8–9% could not be traced back to any still existing star clusters or associations, but are consistent with originating in a cluster which could not be found due to several effects such as incompleteness of their cluster sample (Schilbach & Röser 2008), dissolution of the cluster (Gvaramadze et al. 2012), and having experienced the two-step-ejection process proposed by Pflamm-Altenburg & Kroupa (2010). The majority of the stars which can be traced back to clusters appear to have been ejected when their parental clusters were very young (Schilbach & Röser 2008). This may suggest that dynamical ejections through few-body interactions in the heavy-star-rich cluster core mainly populate the field O star population (Clarke & Pringle 1992).

There have been many studies on dynamically ejected massive stars from young star clusters using numerical integrations. In most cases the studies performed few-body scattering experiments focused on runaways, which are ejected with a velocity higher than ≈30–40 km s⁻¹ (Gualandris et al. 2004; Pflamm-Altenburg & Kroupa 2006; Gvaramadze et al. 2009; Gvaramadze & Gualandris 2011). Full N-body calculations for a particular cluster mass have also been presented (Banerjee et al. 2012; Perets & Šubr 2012). Fujii & Portegies Zwart (2011) performed N-body calculations for very dense star clusters composed of single stars with three different cluster masses ($\geqslant 6000\;{{M}_{\odot }}$) and showed that the clusters with the lowest mass in their models, $6000\;{{M}_{\odot }}$, produce runaways most efficiently. They argued that the one bullying binary formed during core collapse of the cluster is responsible for the production of runaways and that the number of runaways is almost independent of cluster mass. The runaway fraction (number of runaway stars over total number of stars still present) therefore decreases with cluster mass as the more massive clusters have more stars. Note that the lower mass limit for stars in their models is $1\;{{M}_{\odot }}$. Thus, if a cluster is populated with a full range of stellar masses from the hydrogen burning limit upwards, in order to have the same number of massive stars, the cluster would be twice as massive as they claimed. The relation between O star ejection fraction and cluster mass has not been studied in detail to date.

Here we study for the first time how the dynamical ejection fraction of O stars changes with the masses of very young star clusters using direct N-body calculations with a broad range of cluster masses, from low-mass clusters containing one or two O stars to massive clusters containing a few hundred O stars, and including primordial binaries. Such an investigation based on a very large library of young cluster models has not been possible to date because the knowledge on the properties of initial stellar population, their birth configurations, and, last but not least, the computational algorithms enabling such a CPU-intensive and massive numerical challenge have not been in place. Thus, for example, how to initialize initially a mass-segregated cluster in dynamical equilibrium is only known since the work of Šubr et al. (2008) and Baumgardt et al. (2008) and the statistical distribution functions of massive star binaries are becoming observationally constrained only now (e.g., Kiminki & Kobulnicky 2012; Sana et al. 2012, 2013; Kobulnicky et al. 2014; Aldoretta et al. 2015). The N-body calculations, made possible by the relentless code and algorithmic advances by Sverre Aarseth and Seppo Mikkola, used in this study are briefly described in Section 2. The results are presented in Sections 3–5, and the discussion and summary follow in Sections 6 and 7, respectively.

The data used in this study are based on the theoretical young star cluster library of model sequences over cluster mass computed by Oh & Kroupa (2012) with direct N-body calculations using the NBODY6¹ code (Aarseth 2003) and various initial conditions. Among the model sequences in the library, we adopt four sequences here and, further, extend them to higher cluster masses. Additionally we perform calculations for three more sets of initial conditions for this study. Table 1 lists properties of the model sequences studied in this paper. The initial setup of the models is described below.

Table 1.  List of N-body Model Sequences Studied Here

Name	${{r}_{h}}$ (pc)	IMS	IBF	IPD	Pairing method
MS3OP	0.3	Yes	1	Kroupa95	OP
NMS3OP	0.3	No	1	Kroupa95	OP
MS3OP_SP	0.3	Yes	1	Sana et al.12	OP
MS3UQ_SP	0.3	Yes	1	Sana et al.12	Uniform q dist.
MS3S	0.3	Yes	0	⋯	⋯
MS8OP	0.8	Yes	1	Kroupa95	OP
MS1OP	0.1	Yes	1	Kroupa95	OP

Note. Names of model sequences and initial half-mass radii (${{r}_{h}}$) are listed in columns 1 and 2, respectively. Columns 3–6 present initial mass segregation (IMS), initial binary fraction (IBF), and initial period distribution (IPD) and pairing method used for massive binaries (primary mass $\;\geqslant 5\;{{M}_{\odot }}$).

Download table as:  ASCIITypeset image

The four model sequences that we adopt from the library for this study are composed of binary-rich clusters with initial mass segregation (MS3OP) and without initial mass segregation (NMS3OP), single-star clusters (MS3S) with an initial half-mass radius, ${{r}_{h}}$, of 0.3 pc, and initially mass-segregated, binary-rich clusters with an initial ${{r}_{h}}$ of 0.8 pc (MS8OP).

The cluster masses, ${{M}_{{\rm ecl}}}$, in the library range from 10 to ${{10}^{3.5}}\;{{M}_{\odot }}$ (to ${{10}^{4}}\;{{M}_{\odot }}$ for the MS3OP model sequence). For extending the range of cluster masses, we additionally carried out calculations for clusters with ${{M}_{{\rm ecl}}}\approx {{10}^{4}}$ and ${{10}^{4.5}}\;{{M}_{\odot }}$ (${{M}_{{\rm ecl}}}\approx {{10}^{4.5}}\;{{M}_{\odot }}$ for the MS3OP model sequence). Individual stellar masses are randomly drawn from the canonical initial mass function (IMF), which is a two-part power-law (Kroupa 2001; Kroupa et al. 2013),

$$\begin{eqnarray}\xi (m)=k\left\{ \begin{array}{ccccccccccccccc} {{\left( \frac{m}{0.08} \right)}^{-1.3}}, & 0.08\leqslant m/\;{{M}_{\odot }}\lt 0.50, \\ {{\left( \frac{0.5}{0.08} \right)}^{-1.3}}{{\left( \frac{m}{0.5} \right)}^{-2.3}}, & 0.50\leqslant m/\;{{M}_{\odot }}\leqslant {{m}_{{\rm max} }}. \\ \end{array} \right.\end{eqnarray} \tag{ 1 }$$

More details of our stellar mass sampling can be found in Oh & Kroupa (2012). The mass of the most massive star in the cluster, ${{m}_{{\rm max} }}$, is chosen from the maximum-stellar-mass–cluster-mass relation (Weidner & Kroupa 2004; Weidner et al. 2010, 2013a, 2014). Thus only clusters with ${{M}_{{\rm ecl}}}\geqslant {{10}^{2.5}}\;{{M}_{\odot }}$ initially have O-type stars in our models. Pflamm-Altenburg et al. (2007) provide a fitting formula for ${{m}_{{\rm max} }}({{M}_{{\rm ecl}}})$. This procedure was adopted rather than optimal sampling introduced in Kroupa et al. (2013) because the theoretical young cluster library of Oh & Kroupa (2012) did not have this sampling method available.

Binary-rich models have an initial binary fraction of 100%, i.e., all stars are in a binary system at t = 0 Myr. This is motivated by the angular momentum problem of star formation since the angular momentum of a collapsing cloud core can be distributed efficiently into two stars, and also by the strong empirical evidence that star formation universally prefers the binary mode (Goodwin & Kroupa 2005; Duchêne et al. 2007; Goodwin et al. 2007; Marks & Kroupa 2012; Sana et al. 2014; Leigh et al. 2015). The Kroupa (1995b) period distribution function with minimum and maximum periods of 10 and ${{10}^{8.43}}$ days, respectively, is adopted for all binaries in the library. The Kroupa period distribution,

$$\begin{eqnarray}&&{{f}_{P}}({{{\rm log} }_{10}}P)=2.5\frac{{{{\rm log} }_{10}}\;P-1}{45+{{({{{\rm log} }_{10}}\;P-1)}^{2}}},\end{eqnarray} \tag{ 2 }$$

where the period (P) is in days, has been derived iteratively using binary-star data from the Galactic field population and from star-forming regions. With this distribution, about 2.75% of all O star binaries have a period shorter than 100 days (Figure 1). Although it was derived from binaries with primary star mass $\leqslant 1\;{{M}_{\odot }}$, the same distribution is used for more massive binaries in the library since no well constrained period distribution of massive binary population was available when the N-body calculations were begun.

Zoom In Zoom Out Reset image size

However, recently, Sana et al. (2012) derived an intrinsic period distribution of O stars of the form ${{f}_{P}}({{{\rm log} }_{10}}P)\propto {{({{{\rm log} }_{10}}\;P)}^{-0.55}}$, over the period range of $0.15\lt {{{\rm log} }_{10}}(P/{\rm days})\lt 3.5$, from spectroscopic observations of O star populations in nearby young open star clusters. The Sana et al. (2012) period distribution leads to $\approx 49$% of O stars being binaries with a period shorter than 100 days (Figure 1). Here we additionally perform a set of N-body integrations with massive binaries ($m\geqslant 5\;{{M}_{\odot }}$) having the Sana et al. (2012) period distribution, MS3OP_SP, but with the maximum period extended to the value ${{{\rm log} }_{10}}(P/{\rm days})\approx 6.7$ at which the cumulative binary fraction becomes unity. We thus adopt the following period distribution function for the massive binaries,

$$\begin{eqnarray}&&{{f}_{P}}({{{\rm log} }_{10}}P)=0.23\times {{({{{\rm log} }_{10}}\;P)}^{-0.55}},\end{eqnarray} \tag{ 3 }$$

for which $\int _{0.15}^{6.7}{{f}_{P}}({{{\rm log} }_{10}}P)\;d{{{\rm log} }_{10}}\;P=1$. The other initial conditions are the same as those of MS3OP. See also Kiminki & Kobulnicky (2012), Kobulnicky et al. (2014) and Aldoretta et al. (2015) for slightly different period distributions of O star binaries derived from observations.

Stars with a mass $m\lt 5\;{{M}_{\odot }}$ are randomly paired after choosing the two masses randomly from the IMF. Again, this mass-ratio distribution is derived from the Galactic field population and from star-forming regions (Kroupa 1995a, 1995b). Stars with $m\geqslant 5\;{{M}_{\odot }}$ are paired to imprint observations that the massive binaries favor massive companions (Pinsonneault & Stanek 2006; Kobulnicky & Fryer 2007; Kiminki & Kobulnicky 2012; Sana et al. 2012). To achieve this we first order all stars more massive than $5\;{{M}_{\odot }}$ in an array of decreasing mass. The most massive star is then paired with the next massive star, the third most massive star is paired with the fourth most massive star, and so on. We refer to this procedure as ordered pairing (OP). While with this method it is easy to have a massive primary being paired with a massive companion while preserving the IMF, it gives a significant bias toward high mass ratios ($q\approx 1$ where $q={{m}_{2}}/{{m}_{1}}$ and ${{m}_{1}}\geqslant {{m}_{2}}$) which is not supported by recent observations (Kiminki & Kobulnicky 2012; Sana et al. 2012; Kobulnicky et al. 2014). The recent observations show the mass-ratio distribution of O-type binaries to be rather uniform,

$$\begin{eqnarray*}&&f(q)\propto {{q}^{\eta }}\end{eqnarray*}$$

with $-0.1\lesssim \eta \lesssim 0.1$ and $0.1\leqslant q\leqslant 1.0$. We add a model sequence (MS3UQ_SP) with a uniform mass-ratio distribution, $\eta =0$, for massive binaries (${{m}_{1}}\geqslant 5\;{{M}_{\odot }}$), being otherwise the same as the MS3OP_SP model sequence. First, a mass-ratio value for a primary is derived from a uniform probability function and then the secondary is chosen from the rest of the members which gives the closest q value for the derived one. This is important to preserve the IMF: it is important to first generate a full list of stars whose combined mass corresponds to the mass of the cluster and only then to create the binaries using this list only. The thermal distribution is used for the initial eccentricity distribution, i.e., ${{f}_{e}}(e)=2e$, where e is the eccentricity of a binary (see Kroupa 2008).

Each model sequence has one initial half-mass radius regardless of the cluster mass since by observation there is no significant relation between the size and mass of a cluster (Bastian et al. 2005; Scheepmaker et al. 2007) or at most only a very weak one (Larsen 2004; Marks & Kroupa 2012). The cluster-mass–initial half-mass radius relation from Marks & Kroupa (2012),

$$\begin{eqnarray*}&&({{r}_{h}}/{\rm pc})=0.1\times {{({{M}_{{\rm ecl}}}/\;{{M}_{\odot }})}^{0.13}},\end{eqnarray*}$$

gives initial half-mass radii of $\approx 0.21$–0.38 pc for ${{M}_{{\rm ecl}}}={{10}^{2.5}}$–${{10}^{4.5}}\;{{M}_{\odot }}$. Our main, initial half-mass radius (0.3 pc) lies close to the median value and is well consistent with the typical initial half-mass radius (0.1–0.3 pc) inferred for the Milky Way clusters (Marks & Kroupa 2011). Thus, our results may be able to represent reality. Here we include two more model sequences with different half-mass radii of 0.1 and 0.8 pc, MS1OP and MS8OP, respectively. This allows us to estimate how the ejection fraction changes with cluster mass when the cluster radius varies with cluster mass. While the MS8OP model sequence is adopted from the library of Oh & Kroupa (2012), the MS1OP model sequence is computed for this work. From the choices of half-mass radii, the initial central densities of our models range from $\approx 300$ to $3.3\times {{10}^{4}}\;{{M}_{\odot }}\;{\rm p}{{{\rm c}}^{-3}}$ for ${{r}_{h}}=0.8$ pc clusters, from $6.2\times {{10}^{3}}$ to $6.2\times {{10}^{5}}\;{{M}_{\odot }}\;{\rm p}{{{\rm c}}^{-3}}$ for ${{r}_{h}}=0.3$ pc clusters, and from $\approx 1.7\times {{10}^{5}}\;{{M}_{\odot }}\;{\rm p}{{{\rm c}}^{-3}}$ to $\approx 1.7\times {{10}^{7}}\;{{M}_{\odot }}\;{\rm p}{{{\rm c}}^{-3}}$ for ${{r}_{h}}=0.1$ pc clusters.

For all models, initial positions and velocities of single stars or of centers-of-mass of binaries in the clusters are generated following the Plummer density profile under the assumption that the clusters are in virial equilibrium. Initially mass-segregated clusters are produced being fully mass-segregated with the method developed in Baumgardt et al. (2008) in which more massive systems are more bound to the cluster.

The N-body code regularises the motions of stars in compact configurations and incorporates a near-neighbor force evaluation scheme to guarantee the high accuracy of the solutions of the equations of motions. Stellar evolution is taken into account by using the stellar evolution library (Hurley et al. 2000, 2002) implemented in the code. All clusters are evolved up to 3 Myr, i.e., to a time before the first supernova occurs, in order to study dynamical ejections only. Moreover the dynamical ejections of O stars are most efficient at an early age of the cluster (S. Oh et al. 2015, in preparation and see Section 5), thus computations to an older cluster age would not change our conclusion on the dynamical ejections.

A standard solar-neighborhood tidal field is adopted. The clusters are highly tidally underfilling such that the tidal field of the hosting galaxy plays no significant role for the results presented here. Our cluster models are initially gas-free, i.e., without a background gas-potential, and thus the gas-expulsion phase is not included here. While 100 realizations with different random seed numbers are performed for the clusters with ${{M}_{{\rm ecl}}}\leqslant {{10}^{3.5}}\;{{M}_{\odot }}$, 10 and 4 realizations are carried out for ${{M}_{{\rm ecl}}}={{10}^{4}}$ and ${{10}^{4.5}}\;{{M}_{\odot }}$ clusters, respectively, as computational costs increase with the number of stars, i.e., cluster mass, and stochastic effects are reduced with a larger number of stars.

Throughout this paper, an O star refers to a star with a mass of $\geqslant 17.5\;{{M}_{\odot }}$ at 3 Myr and an O star system is either a single O star or a binary system having at least one O star component. We regard a system as being ejected when its distance from the cluster center is greater than $3\times {{r}_{h}}$ of the cluster at 3 Myr, and its velocity is larger than the escape velocity at that radius. Such a short distance criterion for ejections is an appropriate physical condition to quantify the real ejected O star fraction, in order to account for the O stars which have been dynamically ejected but have not travelled far from the cluster center by 3 Myr. O stars can be ejected with a velocity of a few km s⁻¹ from the lower-mass clusters, e.g., the central escape velocity of a cluster with ${{M}_{{\rm ecl}}}=1000\;{{M}_{\odot }}$ and ${{r}_{h}}=0.3$ pc is $\approx 6$ km s⁻¹.

From the computational data the ejection fraction of O star systems, $^{{\rm com}}{{f}_{{\rm ej},{\rm O}}}$, is defined as the number ratio of the ejected O star systems, $^{{\rm com}}{{N}_{{\rm ej},{\rm O}}}$, to all O star systems that are still present in a calculation, $^{{\rm com}}{{N}_{{\rm O}}}$, i.e.,

$$\begin{eqnarray*}&&^{{\rm com}}{{f}_{{\rm ej},{\rm O}}}=\frac{^{{\rm com}}{{N}_{{\rm ej},{\rm O}}}}{^{{\rm com}}{{N}_{{\rm O}}}}.\end{eqnarray*}$$

Only clusters with ${{M}_{{\rm ecl}}}\geqslant {{10}^{2.5}}\;{{M}_{\odot }}$ are shown in the tables because clusters with ${{M}_{{\rm ecl}}}\leqslant 100\;{{M}_{\odot }}$ do not have an O star in our models. Since values of $^{{\rm com}}{{f}_{{\rm ej},{\rm O}}}$ from individual clusters are diverse, especially for low-mass clusters, we mainly deal with the averaged value of $^{{\rm com}}{{f}_{{\rm ej},{\rm O}}}$, ${{\langle }^{{\rm com}}}{{f}_{{\rm ej},{\rm O}}}\rangle$, for the same cluster mass in a certain set of initial conditions.

The average ejection fractions of O star systems for all models are shown in Figure 2. Table A1 lists averaged properties of the models, such as half-mass radius, number of O star systems, number of ejected O star systems, O star ejection fraction, and O star binary fraction in the clusters and among the ejected O star systems, at 3 Myr. All models show a trend that ${{\langle }^{{\rm com}}}{{f}_{{\rm ej},{\rm O}}}\rangle$ increases with cluster mass up to ${{10}^{3.5}}\;{{M}_{\odot }}$ and then decreases for a higher cluster mass. The average values of the ejection fraction with the same cluster mass vary from model to model. The peak fractions are $\approx 25$% for the binary-rich cluster models with ${{r}_{h}}=0.3$ pc with a weak dependency on initial mass segregation or on the initial period distribution of massive binaries, although the ${{\langle }^{{\rm com}}}{{f}_{{\rm ej},{\rm O}}}\rangle$ values of the other cluster masses change with the initial conditions. Note that individual clusters in the mass range ${{10}^{2.5}}$–${{10}^{3.5}}\;{{M}_{\odot }}$ may eject all their O star systems!

Zoom In Zoom Out Reset image size

The ejection fractions from the initially not mass-segregated models (NMS3OP) exhibit similar results to the MS3OP models. The difference between the two models is only seen at the lowest cluster mass, ${{10}^{2.5}}\;{{M}_{\odot }}$. This can be understood because dynamical mass segregation is inefficient for the lowest mass clusters due to their low density. Initial mass segregation therefore does not have a significant impact on massive-star ejections.

In the models using the Sana et al. (2012) period distribution (MS3OP_SP and MS3UQ_SP), the same trend of the O star ejection fraction is present as above, but the ${{\langle }^{{\rm com}}}{{f}_{{\rm ej},{\rm O}}}\rangle$ values show a much broader peak than those in the model using the Kroupa (1995b) period distribution (MS3OP). The ${{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle$ values of the three models are similar (see Table A1). However, for massive clusters (${{M}_{{\rm ecl}}}\gt {{10}^{3.5}}\;{{M}_{\odot }}$), MS3OP_SP clusters not only have a smaller number of O star systems near 3 Myr due to a higher fraction of binaries and mergers, but they also eject more O star systems. The latter suggests that short period massive binaries are dynamically active and boost the ejection of O stars, especially for massive clusters.

The numbers of O star systems are slightly higher in the MS3UQ_SP sequence than in the MS3OP_SP sequence (see Table A1) because O star binaries in the MS3OP_SP models mostly have an O star companion while companions of O star binaries in the MS3UQ_SP models are chosen from a broader range of stellar masses. For example, when the binary fraction is unity and both sequences have the same number of individual O stars, for the MS3OP_SP sequence the number of O star systems would be almost half of the number of individual O stars while it would be larger than half in the case of the MS3UQ_SP sequence. However, the difference in the mass-ratio distributions of the two model sequences does not seem to affect the trend of the O star ejection fraction as a function of cluster mass, having a peak at ${{M}_{{\rm ecl}}}={{10}^{3.5}}\;{{M}_{\odot }}$. The difference between the results of the MS3OP_SP and MS3UQ_SP models is only marginal. Therefore, massive star ejections are not sensitive to the form of the mass-ratio distribution as long as they are paired with other massive stars, i.e., $q\gt 0.1$.

For single-star clusters (MS3S), we find the same cluster mass at which the peak of the O star ejection fraction occurs although the ejection fractions are about half those of the binary-rich model sequences. The peak appears at the same cluster mass in the sequences with different initial cluster sizes. In the case of MS8OP (${{r}_{h}}=0.8$ pc), the ejection fraction is much smaller than in other model sequences with ${{r}_{h}}=0.3$ pc. But ${{\langle }^{{\rm com}}}{{f}_{{\rm ej},{\rm O}}}\rangle$ reaches $\approx 50$% in the MS1OP model sequence (initial ${{r}_{h}}=0.1$ pc) at the cluster mass of the peak. The dependency of the ejection fraction on cluster mass is thus independent of the initial cluster size, but more compact cluster do have higher ejection fractions of massive stars.

The peak near ${{10}^{3.5}}\;{{M}_{\odot }}$ of the ejection fraction can be understood as a consequence of the growth of the number of O stars with cluster mass surpassing the number of ejected O stars as the cluster mass increases and the potential smoothens. Figure 3 shows the average number of ejected O stars at each cluster mass.

Zoom In Zoom Out Reset image size

This implies that clusters a few Myr old with a mass at which the O star ejection fraction peaks would present O star populations most deviating from the number of O stars given by the canonical IMF. For example the estimated stellar mass in the inner 2 pc of the Orion Nebula Cluster (ONC) is $\approx 1800\;{{M}_{\odot }}$ (Hillenbrand & Hartmann 1998) and probably $\approx 2700\;{{M}_{\odot }}$ taking account of a probable binary fraction of 50% (Kroupa 2000). With these cluster masses, about 6–10 stars are expected to have a mass larger than $18\;{{M}_{\odot }}$ from the canonical IMF. However, only three O-type (or $m\geqslant 17.5\;{{M}_{\odot }}$) stars are found within $\approx 2$ pc of the Trapezium stars (Hillenbrand 1997). This discrepancy may be due to the stochastic sampling of the IMF. But the probability for this is less than 5% when stellar masses are randomly drawn from the IMF with a stellar upper mass limit of $150\;{{M}_{\odot }}$ (but about 10% with ${{m}_{{\rm max} }}=50\;{{M}_{\odot }}$) such that it is more likely that the ONC has ejected a significant fraction of its O star content. Furthermore, the initial mass of the cluster may have been higher than observed at the present day as the cluster could have lost a significant fraction of its members during the residual gas expulsion phase (Kroupa et al. 2001). Thus the ONC formed with a mass close to the peak mass we found in this study and could have lost a significant fraction of its massive stars through the dynamical ejection process. This may explain why the ONC's massive star population deviates from the canonical IMF (Pflamm-Altenburg & Kroupa 2006). Examples of such events would be two O-type single runaway stars, AE Aur and μ Col, and a massive binary system ι Ori which may have been ejected from the ONC through a binary–binary interaction (Gies & Bolton 1986; Hoogerwerf et al. 2001; Gualandris et al. 2004). Future proper motion surveys, such as with the Gaia - mission, may provide further constraints on this question, because the putative ejected O stars from the ONC would have to become evident.

Very massive clusters may have small ejection fractions of O stars but the upper end of their stellar mass range could lack stars because the ejection fraction of stars more massive than $100\;{{M}_{\odot }}$ can be from $\approx 10$% to $\approx 100$% (Banerjee et al. 2012). Furthermore, they can eject not only very massive single stars but also very massive binaries (Oh et al. 2014). This bias through the most massive stars being ejected from star-burst clusters and the observed IMF of star-burst clusters being close to the canonical value (e.g., the R136 cluster in the Large Magellanic Cloud, Massey & Hunter 1998) may imply the IMF to be top-heavy in star-burst clusters (Banerjee & Kroupa 2012). There is some evidence for a possibly systematically varying IMF above a few ${{M}_{\odot }}$: a number of independent arguments have shown the IMF to become increasingly top-heavy with increasing star-formation rate density on a pc scale (Marks et al. 2012), with a dependency on the metallicity. This metallicity dependence is such that metal-rich environments tend to produce a steeper (i.e., top-light) IMF slope. Evidence for this may have emerged in M31 data (Weisz et al. 2015). The models calculated here are, however, in the invariant regime.

A theoretically expected relation between the ejection fraction and the cluster mass can be found. First, the number of O stars in a cluster, ${{N}_{{\rm O},{\rm IMF}}}$, which is a function of ${{m}_{{\rm max} }}({{M}_{{\rm ecl}}})$ and ${{M}_{{\rm ecl}}}$, is

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{N}_{{\rm O},{\rm IMF}}} & = & \int _{17.5\,{{M}_{\odot }}}^{{{m}_{{\rm max} }}({{M}_{{\rm ecl}}})}\xi (m)dm, \\ {} & \approx & \;\frac{3.390\times {{10}^{-4}}-0.014m_{{\rm max} }^{-1.3}}{0.100-0.062m_{{\rm max} }^{-0.3}}{{M}_{{\rm ecl}}}, \\ \end{array}\end{eqnarray} \tag{ 4 }$$

by using the IMF (Equation (1)). The normalization constant, k, in Equation (1) is calculated from ${{M}_{{\rm ecl}}}=\int _{0.08\,{{M}_{\odot }}}^{{{m}_{{\rm max} }}}m\xi (m)dm.$ The top panel of Figure 4 presents ${{N}_{{\rm O},{\rm IMF}}}$ as a function of ${{M}_{{\rm ecl}}}$ from Equation (4). Since we adopt the ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation, Equation (4) is not a simple linear function of cluster mass especially for low-mass clusters in which ${{m}_{{\rm max} }}$ is highly dependent on the cluster mass. ${{N}_{{\rm O},{\rm IMF}}}$, however, becomes almost linearly proportional to ${{M}_{{\rm ecl}}}$ for massive clusters where ${{m}_{{\rm max} }}$ is saturated (see the upper panel of Figure 4). Note that the solid line in the upper panel of Figure 4 and the numbers in Table A1 are slightly different because Equation (4) gives the number of all individual O stars from the IMF, whereas ${{\langle }^{{\rm com}}}{{N}_{{\rm O}}}\rangle$ listed in Table A1 is the number of O star systems at 3 Myr. The difference between the number of O stars from the IMF, ${{N}_{{\rm O},{\rm IMF}}},$ and from our N-body calculations, ${{\langle }^{{\rm com}}}{{N}_{{\rm O}}}\rangle$, at the cluster mass of ${{10}^{2.5}}\;{{M}_{\odot }}$ in Figure 4 arises due to our choice of forcing a cluster to have at least one ${{m}_{{\rm max} }}$ star initially. This would not have been required if optimal sampling (Kroupa et al. 2013) had been available when the theoretical young star cluster library was computed pre-2013.

Zoom In Zoom Out Reset image size

Our stellar masses are randomly drawn from the IMF with the upper mass limit given by the ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation. Therefore, our sampling would produce a slightly lower maximum stellar mass in a cluster but a similar number of O stars compared to what random sampling without the ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation would give. The gray line in the upper panel of Figure 4 is the expected number of O stars from the IMF with a universal upper stellar mass limit of $150\;{{M}_{\odot }}$ for all cluster masses. The figure indicates that a higher number of O stars is expected in the IMF with a universal and constant ${{m}_{{\rm max} }}=150\;{{M}_{\odot }}$ than with the ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation at lower cluster mass (${{M}_{{\rm ecl}}}\;\lesssim$ a few $1000\;{{M}_{\odot }}$, the black solid line in the same figure). But the values of ${{\langle }^{{\rm com}}}{{N}_{{\rm O}}}\rangle$ in our N-body models are close to the gray line because of adding one ${{m}_{{\rm max} }}({{M}_{{\rm ecl}}})$ star in a cluster of mass ${{M}_{{\rm ecl}}}$ by default. So whether adopting the relation or not would have little effect on the number of O stars in our models. The proper comparison for outcomes of the two assumptions (${{m}_{{\rm max} }}=150\;{{M}_{\odot }}$ or the ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation) can be done by performing a large set of Monte Carlo calculations. However, that is beyond the objectives of this study. In any case, the ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation is observationally well supported for young star clusters and star-forming regions in the Milky Way given the most recent data (Kirk & Myers 2011; Weidner et al. 2013a).

Fitting the number of O star systems present, regardless of whether the systems are in the star cluster or ejected, in N-body calculations at 3 Myr to a simple power law gives,

$$\begin{eqnarray}{{\langle }^{{\rm com}}}{{N}_{{\rm O}}}\rangle \propto \left\{ \begin{array}{ccccccccccccccc} M_{{\rm ecl}}^{0.97\pm 0.01}, & {\rm MS}3{\rm OP}, \\ M_{{\rm ecl}}^{0.96\pm 0.01}, & {\rm NMS}3{\rm OP}, \\ M_{{\rm ecl}}^{0.96\pm 0.01}, & {\rm MS}3{\rm OP}\_{\rm SP}, \\ M_{{\rm ecl}}^{0.96\pm 0.01}, & {\rm MS}3{\rm UQ}\_{\rm SP}, \\ M_{{\rm ecl}}^{1.01\pm 0.01}, & {\rm MS}3{\rm S}, \\ M_{{\rm ecl}}^{0.95\pm 0.01}, & {\rm MS}8{\rm OP}, \\ M_{{\rm ecl}}^{1.00\pm 0.01}, & {\rm MS}1{\rm OP}. \\ \end{array} \right.\end{eqnarray} \tag{ 5 }$$

The numbers of O star systems in all model sequences are more or less linearly proportional to the cluster mass as expected from the IMF despite various binary fractions.

Next, the relation between the number of ejected O stars and cluster mass is needed. According to Verbunt (2003), the binary–single star scattering encounter rate, Γ, is

$$\begin{eqnarray}&&{\Gamma }\propto \frac{\rho _{0}^{2}r_{c}^{3}}{{{{\boldsymbol{ \sigma }} }_{{\rm cl}}}}a,\end{eqnarray} \tag{ 6 }$$

where ${{\rho }_{0}}$ is the central mass density, ${{\sigma }_{{\rm cl}}}$ the velocity dispersion of the cluster, ${{r}_{c}}$ the core radius, and a the semimajor axis of the binary. The semi-major axis can be replaced with the system mass, ${{m}_{{\rm sys}}}$, and the orbital velocity, ${{v}_{{\rm orb}}}$, of the binary by using Equation (8.142) in Kroupa (2008) as,

$$\begin{eqnarray}&&a\propto {{m}_{{\rm sys}}}v_{{\rm orb}}^{-2}.\end{eqnarray} \tag{ 7 }$$

Considering that binaries with a circular velocity similar to the velocity dispersion of the cluster are dynamically the most active (i.e., ${{v}_{{\rm orb}}}={{\sigma }_{{\rm cl}}}$) and ${{\sigma }_{{\rm cl}}}\propto \sqrt{{{\rho }_{0}}}{{r}_{c}}$, Equation (6) becomes

$$\begin{eqnarray}&&{\Gamma }\propto \rho _{0}^{0.5}{{m}_{{\rm sys}}}\propto M_{{\rm ecl}}^{0.5}{{m}_{{\rm sys}}},\end{eqnarray} \tag{ 8 }$$

since all our cluster models in a given sequence have the same initial half-mass radius so the initial central density is only proportional to the cluster mass. A binary ejecting an O star after a strong interaction with a single star is likely a binary consisting of O stars since the least massive star is generally ejected after a close encounter between a binary and a single star. Therefore, ${{m}_{{\rm sys}}}$ can be approximated to twice the average mass of O stars in a cluster. We assume,

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{m}_{{\rm sys}}} & = & 2\times \langle {{m}_{\geqslant 17.5\,{{M}_{\odot }}}}\rangle =2\frac{\int _{17.5}^{{{m}_{{\rm max} }}}m\xi (m)dm}{\int _{17.5}^{{{m}_{{\rm max} }}}\xi (m)dm} \\ {} & \approx & 8.66\frac{m_{{\rm max} }^{-0.3}-0.42}{m_{{\rm max} }^{-1.3}-0.02}, \\ \end{array}\end{eqnarray} \tag{ 9 }$$

using the IMF (Equation (1)). Thus ${{m}_{{\rm sys}}}$ is a function of ${{m}_{{\rm max} }}$ so that ${{m}_{{\rm sys}}}={{m}_{{\rm sys}}}({{m}_{{\rm max} }})$. Due to the ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation, ${{m}_{{\rm sys}}}$ is thus a function of the cluster mass (the lower panel of Figure 4).

The number of ejections is proportional to the encounter rate,

$$\begin{eqnarray}&&^{{\rm th}}{{N}_{{\rm ej},{\rm O}}}\propto {\Gamma }\propto {{m}_{{\rm sys}}}({{m}_{{\rm max} }})M_{{\rm ecl}}^{0.5}.\end{eqnarray} \tag{ 10 }$$

The ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation has only a minor effect on the number of ejected O stars for the massive clusters, i.e., the dependence of ${{m}_{{\rm sys}}}$ on ${{M}_{{\rm ecl}}}$ becomes negligible in Equation (10). $^{{\rm th}}{{N}_{{\rm ej},{\rm O}}}$ is thus expected to be proportional to $M_{{\rm ecl}}^{0.5}$ for the massive clusters. We fit ${{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle$ as a power law of cluster mass for clusters with ${{10}^{3.5}}\leqslant {{M}_{{\rm ecl}}}/\;{{M}_{\odot }}\leqslant {{10}^{4.5}}$ (Figure 3). These fits provide for each model sequence,

$$\begin{eqnarray}{{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle \propto \left\{ \begin{array}{ccccccccccccccc} M_{{\rm ecl}}^{0.63\pm 0.066}, & {\rm MS}3{\rm OP}, \\ M_{{\rm ecl}}^{0.67\pm 0.066}, & {\rm NMS}3{\rm OP}, \\ M_{{\rm ecl}}^{0.78\pm 0.034}, & {\rm MS}3{\rm OP}\_{\rm SP}, \\ M_{{\rm ecl}}^{0.76\pm 0.061}, & {\rm MS}3{\rm UQ}\_{\rm SP}, \\ M_{{\rm ecl}}^{0.62\pm 0.061}, & {\rm MS}3{\rm S}, \\ M_{{\rm ecl}}^{0.66\pm 0.100}, & {\rm MS}8{\rm OP}, \\ M_{{\rm ecl}}^{0.87\pm 0.100}, & {\rm MS}1{\rm OP}. \\ \end{array} \right.\end{eqnarray} \tag{ 11 }$$

All model sequences but MS3OP_SP, MS3UQ_SP, and MS1OP exhibit slopes similar to each other within an error, in a range of 0.62–0.67. The power-law exponents we obtain from our N-body models are slightly steeper than 0.5 in Equation (10) but with only a marginal difference. It can be argued that this is due to a broad range of binary parameters for which binaries could have engaged in the ejection of O stars in our models and the complexity of the few-body scattering process. The model sequences MS3OP_SP and MS3UQ_SP show steeper slopes, i.e., the more massive clusters in this model eject O stars more efficiently compared to other models. As the model sequence MS1OP begins with an extreme density the results of the model significantly deviate from the simplified theoretical expectation with simple approximations.

Finally, the theoretically expected ejection fraction of O stars is

$$\begin{eqnarray}&&^{{\rm th}}{{f}_{{\rm ej},{\rm O}}}=\frac{^{{\rm th}}{{N}_{{\rm ej},{\rm O}}}}{{{N}_{{\rm O},{\rm IMF}}}}\propto F({{m}_{{\rm max} }})M_{{\rm ecl}}^{-0.5},\end{eqnarray} \tag{ 12 }$$

where $F({{m}_{{\rm max} }})$ comes from the part of Equations (4) and (10) in which ${{m}_{{\rm max} }}$ plays a role. For massive clusters (${{M}_{{\rm ecl}}}\geqslant {{10}^{3.5}}\;{{M}_{\odot }}$), $F({{m}_{{\rm max} }})$ in Equation (12) becomes very weakly dependent on the cluster mass as in Equation (10), and thus this term can be ignored. Therefore, $^{{\rm th}}{{f}_{{\rm ej},{\rm O}}}$ is expected to be proportional to $M_{{\rm ecl}}^{-0.5}$ in this cluster mass regime. Figure 5 shows the linear fits to the O star ejection fractions as a function of cluster mass for the massive clusters (${{10}^{3.5}}\leqslant {{M}_{{\rm ecl}}}/\;{{M}_{\odot }}\leqslant {{10}^{4.5}}$) in a log–log scale. Within a cluster mass range ${{10}^{3.5}}$–${{10}^{4.5}}\;{{M}_{\odot }}$, the numerical data suggest the relation between O star ejection fraction and cluster mass as follows,

$$\begin{eqnarray}{{\langle }^{{\rm com}}}{{f}_{{\rm ej},{\rm O}}}\rangle \propto \left\{ \begin{array}{ccccccccccccccc} M_{{\rm ecl}}^{-0.40\pm 0.064}, & {\rm MS}3{\rm OP}, \\ M_{{\rm ecl}}^{-0.40\pm 0.055}, & {\rm NMS}3{\rm OP}, \\ M_{{\rm ecl}}^{-0.25\pm 0.029}, & {\rm MS}3{\rm OP}\_{\rm SP}, \\ M_{{\rm ecl}}^{-0.20\pm 0.060}, & {\rm MS}3{\rm UQ}\_{\rm SP}, \\ M_{{\rm ecl}}^{-0.43\pm 0.055}, & {\rm MS}3{\rm S}, \\ M_{{\rm ecl}}^{-0.31\pm 0.094}, & {\rm MS}8{\rm OP}, \\ M_{{\rm ecl}}^{-0.16\pm 0.036}, & {\rm MS}1{\rm OP}. \\ \end{array} \right.\end{eqnarray} \tag{ 13 }$$

The results of most models are reasonably consistent with the relation between the ejection fraction and cluster mass expected from the binary–single star scattering encounter rate (Equation (12)) when applied to massive clusters in which the maximum stellar mass is less sensitive to cluster mass.

Zoom In Zoom Out Reset image size

For lower mass clusters, the number of O stars in a cluster and the ejection rate do not simply follow from the above equations due to the dependency of ${{m}_{{\rm max} }}$ on the cluster mass and the highly collisional nature of the cluster cores which contain but a few O stars.

The ejection fractions from individual clusters are plotted with open circles in Figure 2. For lower mass clusters the O star ejection fraction varies largely from one cluster to another, e.g., the values spread from 0 to 1 for clusters with ${{M}_{{\rm ecl}}}\leqslant {{10}^{3.5}}\;{{M}_{\odot }}$. The spread gets smaller with increasing cluster mass as the stochastic effect diminishes with a larger number of O stars in more massive clusters. For massive clusters (${{M}_{{\rm ecl}}}\geqslant {{10}^{4}}\;{{M}_{\odot }}$) the individual ejection fractions mostly lie close to the average value.

Only for a comparison, we additionally show results using 10 pc as the ejection criterion in Figure 2 (without a constraint on the velocity of the stars) since a larger distance criterion is more appropriate from an observational point of view. The ejection fractions using 10 pc as a distance criterion show a similar shape and the peak at the same cluster mass, only the values are smaller (max. up to $\approx 17$% for clusters with an initial ${{r}_{h}}=0.3$ pc) than the results using our main ejection criteria. The difference becomes smaller or almost disappears for the massive clusters. Massive clusters may eject O stars most efficiently at an earlier time of the evolution due to their short dynamical (i.e., crossing) time-scale, and soon their ejection efficiency would drop as the O star population is reduced in the core of the cluster. Also, in the massive clusters the velocity of the ejected systems is generally high (see Section 5.1) so they quickly move further than 10 pc away from the cluster center. Low-mass clusters have longer dynamical time-scales and they eject O stars moving relatively slowly so that they need more time to reach a distance of 10 pc.

Our N-body results show that ${{{\rm log} }_{10}}\;{{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle$ steeply increases with ${{{\rm log} }_{10}}\;{{M}_{{\rm ecl}}}$ at a lower cluster-mass range (${{M}_{{\rm ecl}}}\lt {{10}^{3.5}}\;{{M}_{\odot }}$) while it can be approximated to a linear function of ${{{\rm log} }_{10}}\;{{M}_{{\rm ecl}}}$ at a higher cluster-mass range. In the previous section, we use data only from the clusters with ${{10}^{3.5}}\;{{M}_{\odot }}\leqslant {{M}_{{\rm ecl}}}\leqslant {{10}^{4.5}}\;{{M}_{\odot }}$ to fit a linear relation between ${{{\rm log} }_{10}}{{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle$ and ${{{\rm log} }_{10}}\;{{M}_{{\rm ecl}}}$ (Figure 3). In order to fit the number of ejected O stars for the entire cluster mass range in which ${{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle$ is $\gt 0$, we use the following functional form,

$$\begin{eqnarray}&&{{{\rm log} }_{10}}\;y=a+b{{{\rm log} }_{10}}\;x-\frac{1}{{{({{{\rm log} }_{10}}\;x-{{{\rm log} }_{10}}\;{{x}_{0}})}^{n}}},\end{eqnarray} \tag{ 14 }$$

where $x\gt {{x}_{0}}$, $y={{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle$, and $x={{M}_{{\rm ecl}}}/\;{{M}_{\odot }}$. In this function, ${{{\rm log} }_{10}}\;y$ converges to $-\infty$ for $x\to {{x}_{0}}$ and approximates a linear function, $a+b{{{\rm log} }_{10}}\;x$, for $x\gg {{x}_{0}}$. This function can thus be fitted to our results. The value x₀ indicates the cluster mass where the y value, i.e., ${{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle$, $\approx 0$. For all models but MS3S and MS8OP, we choose ${{x}_{0}}=240\;{{M}_{\odot }}$ which is the maximum cluster mass at which clusters do not host single O stars initially, as given by the ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation. For the other two models we choose ${{x}_{0}}={{10}^{2.5}}\;{{M}_{\odot }}$, at which mass a cluster contains O stars but ${{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle$ is 0 in the N-body calculations. By performing nonlinear least squares fitting of Equation (14) to our results, we obtain the parameters listed in Table 2 and plot the results in Figure 6.

Zoom In Zoom Out Reset image size

Table 2.  Parameters by Fitting Equation (14) for all Models

Model	a	b	n
MS3OP	−0.56 ± 0.25	0.53 ± 0.07	0.29 ± 0.04
NMS3OP	−0.43 ± 0.25	0.49 ± 0.07	0.42 ± 0.05
MS3OP_SP	−0.91 ± 0.16	0.64 ± 0.04	0.26 ± 0.04
MS3UQ_SP	−1.04 ± 0.25	0.67 ± 0.07	$0.22\pm 0,04$
MS3S	0.67 ± 0.32	0.14 ± 0.08	0.90 ± 0.10
MS8OP	−0.60 ± 0.01	0.37 ± 0.003	0.51 ± 0.03
MS1OP	−1.17 ± 0.16	0.80 ± 0.05	0.16 ± 0.03

Download table as:  ASCIITypeset image

The mass distribution of young star clusters can be described with a simple power law,

$$\begin{eqnarray}&&{{\xi }_{{\rm ecl}}}({{M}_{{\rm ecl}}})={{k}_{{\rm ecl}}}\;M_{{\rm ecl}}^{-\beta },\end{eqnarray} \tag{ 15 }$$

where the normalization constant ${{k}_{{\rm ecl}}}$ is calculated from the total mass of all clusters which form within a certain time. Assuming a global Milky Way SFR of $3\;{{M}_{\odot }}\;{\rm y}{{{\rm r}}^{-1}}$, a minimum cluster mass ${{M}_{{\rm ecl},{\rm min} }}=10\;{{M}_{\odot }}$, and $\beta =2$ for the Milky Way, the maximum cluster mass, ${{M}_{{\rm ecl},{\rm max} }}$, that can form in the Galaxy becomes $\approx 1.9\times {{10}^{5}}\;{{M}_{\odot }}$ (Weidner et al. 2004; Marks & Kroupa 2011, and references therein). The total mass formed in stars over time $\delta t$, M_SCS, is

$$\begin{eqnarray}&&{{M}_{{\rm SCS}}}={\rm SFR}\times \delta t=\int _{{{M}_{{\rm ecl},{\rm min} }}}^{{{M}_{{\rm ecl},{\rm max} }}}{{M}_{{\rm ecl}}}\;{{\xi }_{{\rm ecl}}}({{M}_{{\rm ecl}}})\;d{{M}_{{\rm ecl}}}\end{eqnarray} \tag{ 16 }$$

where $\delta t$ is the cluster-population formation time-scale, about 10 Myr (Weidner et al. 2004). We adopt here the ${{M}_{{\rm ecl},{\rm max} }}$–SFR relation derived by Weidner et al. (2004), i.e., the mass of the most massive forming cluster is dependent on the physical properties of the star-formation environment provided by a self-regulated galaxy. It has been found that the luminosity of the brightest cluster increases with the galaxy-wide SFR (Larsen 2002, and references therein) or with the total number of clusters (Whitmore 2003). The relation is treated to be the result of the size of the sample in several studies (Whitmore 2003; see also Section 2.4 of Portegies Zwart et al. 2010, and references therein). Recently, Pflamm-Altenburg et al. (2013) ruled out that the relation results from the pure size-of-sample effect by studying a radial dependency of maximum star cluster masses in M33. Randriamanakoto et al. (2013) confirm this by studying the ${{M}_{{\rm ecl},{\rm max} }}$ versus SFR relation in their independent observational survey. The authors noted that their result can be explained with purely random sampling but only if the luminosity function of super star clusters at the bright end is very steep (a slope $\gt 2.5$), steeper than usually observed ($\approx 2$).

By combining the fitting function for the number of ejected O stars (Equation (14)) with the cluster mass function (Equation (15)), the number of O stars being ejected to the field contributed by clusters in the mass range ${{M}_{{\rm ecl}}}$ to ${{M}_{{\rm ecl}}}+d{{M}_{{\rm ecl}}}$ is

$$\begin{eqnarray}&&{{d}^{{\rm total}}}{{N}_{{\rm ej},{\rm O}}}({{M}_{{\rm ecl}}})={{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle ({{M}_{{\rm ecl}}})\times {{\xi }_{{\rm ecl}}}({{M}_{{\rm ecl}}})\;d{{M}_{{\rm ecl}}},\end{eqnarray} \tag{ 17 }$$

and is shown in the upper panel of Figure 7. As the number of clusters decreases more steeply (${{\xi }_{{\rm ecl}}}({{M}_{{\rm ecl}}})\propto M_{{\rm ecl}}^{-2}$) than the number of ejected O star systems from a cluster increases with cluster mass (${{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle ({{M}_{{\rm ecl}}})\propto M_{{\rm ecl}}^{\gamma }$ with $0.6\lt \gamma \lt 0.9$, Section 3), the number of O stars contributed to the field by dynamical ejections is a function decreasing with increasing cluster mass. The most prominent contributors to the population of ejected O stars are ≈400 $\;{{M}_{\odot }}$ clusters.

Zoom In Zoom Out Reset image size

In the lower panel of Figure 7, cumulative numbers of ejected O-type systems as a function of cluster mass are shown. The fraction (right Y-axis) is calculated by dividing the cumulative number by the total number of O stars formed during the cluster-population forming time-scale of 10 Myr, which can be deduced from Equations (4) and (15),

$$\begin{eqnarray}&&^{{\rm total}}{{N}_{{\rm O}}}=\int _{{{M}_{{\rm ecl},{\rm min} }}}^{{{M}_{{\rm ecl},{\rm max} }}}{{N}_{{\rm O},{\rm IMF}}}({{M}_{{\rm ecl}}})\;{{\xi }_{{\rm ecl}}}({{M}_{{\rm ecl}}})\;d{{M}_{{\rm ecl}}}.\end{eqnarray} \tag{ 18 }$$

The parameters which we adopt here for the Milky Way provide $^{{\rm total}}{{N}_{{\rm O}}}\approx 66749$ being formed every 10 Myr. Note $^{{\rm total}}{{N}_{{\rm O}}}$ is the total number of individual O stars which have formed within the formation time-scale, while ${{d}^{{\rm total}}}{{N}_{{\rm ej},{\rm O}}}({{M}_{{\rm ecl}}})$ is the number of ejected O star systems. It is difficult to derive the total number of O star systems from the models because a binary fraction evolves due to dynamical disruptions/captures and mergers between stars, the final binary fraction varies with cluster mass and initial conditions, and some O stars have an O star secondary while some do not. The fraction shown in the lower panel of Figure 7 can, thus, be interpreted as being a minimum value since the total number of O star systems would be smaller than $^{{\rm total}}{{N}_{{\rm O}}}$ which counts all individual O stars from the IMF and cluster mass function. Our models thus yield as a lower boundary that binary-rich clusters with an initial ${{r}_{h}}=0.3$ pc can eject ≈13–16% of the total O star systems and even $\approx 38\%$ can be dynamically ejected if clusters form more compactly, e.g., with an initial ${{r}_{h}}=0.1$ pc (Figure 7). About 50% of ejected O stars come from clusters with ${{M}_{{\rm ecl}}}\lesssim 3000$–$5000\;{{M}_{\odot }}$.

The fraction of O stars found in clusters and/or OB associations is 70–76% (de Wit et al. 2005; Eldridge et al. 2011, depending on how to count a binary) yielding the fraction of the field O stars as being 24–30%. Our model is thus reasonably consistent with these data. Additional three processes that are not part of the present model but which increase the apparent number of field O stars are: (i) some of the observed field stars may originate from clusters which contained only a few or one O stars and which dissolved at an early age by rapid gas removal (Gvaramadze et al. 2012). That is, such O stars "lost" their birth cluster through disruptive gas removal (see also Kroupa & Boily 2002; Weidner et al. 2007). (ii) Small-N clusters for which the two-body relaxation time is comparable to the crossing time also dissolve within a few initial crossing times (Moeckel et al. 2012), which is aided by stellar evolution within a few tens of Myr. (iii) Some O stars may be ejected when they are companions to primaries that explode as supernovae at an age older than 3 Myr (Eldridge et al. 2011).

In observations about 6% of O stars are found being runaways which have a velocity higher 30 km s⁻¹ (Eldridge et al. 2011). The runaway fraction among the dynamically ejected O stars generally increases with cluster mass (upper panel of Figure 8), e.g., 0% (${{M}_{{\rm ecl}}}={{10}^{2.5}}\;{{M}_{\odot }}$) to 43% (${{M}_{{\rm ecl}}}={{10}^{4.5}}\;{{M}_{\odot }}$) in the case of the MS3UQ_SP model and 81% for the ${{10}^{5}}\;{{M}_{\odot }}$ model in Banerjee et al. (2012). We deduce the number of runaways first by fitting a function for the runaway fraction as a function of cluster mass with the same functional form of the number of ejected O stars (Equation (14), but y being the runaway fraction) and then by multiplying the fitting function by the number of ejected O stars (Equation (17), upper panel of Figure 7). By dividing the total number of runaways by the total number of ejected stars (lower panel of Figure 7), we obtain the runaway fraction among the ejected O stars. Our most realistic model sequence, MS3UQ_SP, predicts about 30% of the ejected O stars, i.e., about 5% of all O star systems, to be runaways. Eldridge et al. (2011) studied runaways assuming the binary supernova scenario, in which the companion star obtains a high velocity when one of the binary components undergoes a supernova explosion. They predicted about 0.5%² of O stars being runaways when only this process is accounted for. The value combining ours and Eldridge et al. (2011) is thus well consistent with the observed one.

Zoom In Zoom Out Reset image size

The upper panel of Figure 7 demonstrates that the presence or absence of initial mass segregation leads to differences in the numbers of ejected O stars only for lower mass clusters (${{M}_{{\rm ecl}}}\lesssim {{10}^{3}}\;{{M}_{\odot }}$) with the constant half-mass radius we assume in this study. If larger ${{r}_{h}}$ are assumed the cluster mass at which the difference vanishes would move up to a larger cluster mass since the mass-segregation time-scale gets longer.

In this section, we study how properties such as velocities, distances from the cluster center, and system masses of the ejected O star systems vary with cluster mass. We remind the reader that our analysis is restricted to systems younger than 3 Myr.

5.1. Kinematics of the Ejected O Star Systems

Cumulative distributions of velocity, distance from the cluster center, and system mass of all ejected O star systems for each cluster mass and for each set of cluster models are shown in Figure 9. Overall, the systems ejected from more massive clusters have velocity distributions skewed toward higher velocity. This can be understood in two ways. First, the velocity dispersion of stars in a cluster increases with cluster mass under the assumption of a constant cluster size (in general, velocity dispersion is a function of both size and mass of a cluster), i.e., on average stars in more massive clusters have higher velocities. Second, a star needs a velocity greater than the escape velocity, which also increases with cluster mass, in order to leave its parental cluster. For the MS8OP model sequence (with ${{r}_{h}}=0.8$ pc), the trend is not as pronounced as in other models since the ejection fractions are small and so are the total number of ejected O star systems, e.g., only in total 11 O stars are ejected out of 4 runs in the ${{10}^{4.5}}\;{{M}_{\odot }}$ cluster models. Also in the case of the ${{10}^{2.5}}\;{{M}_{\odot }}$ clusters in the NMS3OP sequence (${{r}_{h}}=0.3$ pc) only three O star systems are ejected out of 100 runs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The highest velocity of an ejected system generally increases with cluster mass, although this is not always the case, e.g, in the MS3OP model sequence, the system moving the fastest (a single star of mass $46\;{{M}_{\odot }}$ with a velocity of $\approx 200$ km s⁻¹) is ejected from a ${{10}^{3.5}}\;{{M}_{\odot }}$ cluster. The velocity distribution is dependent on the cluster model, especially on the initial half-mass radius of the cluster. The MS1OP models have velocity distributions biased toward higher velocities compared to the larger-sized cluster models (Figure 9). The highest velocity of an ejected system (a single star of mass $34\;{{M}_{\odot }}$ from an ${{M}_{{\rm ecl}}}={{10}^{4.5}}\;{{M}_{\odot }}$ cluster) is $\approx 300$ km s⁻¹ in the MS1OP model sequence.

As the distance distribution of the ejected systems is directly related to the velocity distribution, their shapes are similar. The slight difference between them results from the time of ejections occurring at systematically different times in different models. Owing to the systematic shift to larger ejection velocities, the ejected systems travel on average further away from the clusters in the more massive clusters. The distances of the ejected systems from ${{M}_{{\rm ecl}}}\geqslant {{10}^{4}}\;{{M}_{\odot }}$ clusters range from a few pc to a few 100 pc. The most extreme system (a single star with a mass of $34\;{{M}_{\odot }}$) ejected from an MS1OP model traveled up to $\approx 750$ pc away from its birth cluster by 3 Myr.

Cumulative mass distributions of the ejected O stars are plotted in the right column of Figure 9. In the case of a binary system, masses of both components are counted separately. As binary fractions of the ejected systems are low for most of the models, the cumulative distributions of system masses would be similar to those of individual stellar masses. A few models suggest that the more massive clusters may have a higher fraction of more-massive ejected stars. However the trend is not strongly evident. For example, in the case of the MS1OP models, the mass distributions of the ejected O stars look similar to each other for different cluster masses. But only clusters with ${{M}_{{\rm ecl}}}\geqslant {{10}^{3.5}}\;{{M}_{\odot }}$ eject O stars heavier than $100\;{{M}_{\odot }}$ as the lower mass clusters do not have such a massive star initially, e.g., for ${{M}_{{\rm ecl}}}=1000\;{{M}_{\odot }}$, ${{m}_{{\rm max} }}\approx 44\;{{M}_{\odot }}$. About 10% of the ejected O stars from clusters with ${{M}_{{\rm ecl}}}\gtrsim {{10}^{3.5}}\;{{M}_{\odot }}$ in the MS3OP_SP and MS3UQ_SP sequences are more massive than $100\;{{M}_{\odot }}$, while this is the case for $\lesssim 5$% from the other models. This may result from a higher number of mergers in the MS3OP_SP and MS3UQ_SP sequences due to the high fraction of close, eccentric binaries in comparison to the other models.

Thin gray lines in the right column of Figure 9 indicate cumulative mass distributions of O stars from the canonical IMF with lower and upper mass limits of, respectively, $17.5\;{{M}_{\odot }}$ and the mass of the most massive star among the ejected O star for each cluster-mass set. It follows that most models show that ejected stars have a shallower slope of the mass function than the canonical value (index of 2.3, see Equation (1)) as also shown in Fujii & Portegies Zwart (2011), resembling the mass function in the cluster core.

The normalized cumulative number of the ejected O star systems as a function of time for a few models is shown in Figure 10. In the case of the MS3OP_SP model (${{M}_{{\rm ecl}}}={{10}^{3.5}}\;{{M}_{\odot }}$), almost 70% of the ejected stars travel beyond three times the half-mass radius before 1.5 Myr. Initially not mass-segregated cluster models (e.g., the ${{10}^{3.5}}\;{{M}_{\odot }}$ clusters of the NMS3OP model sequence in Figure 10) show a delay of the ejections as the model requires time for massive stars to migrate to the cluster center to interact with other massive stars. But about half of the ejected stars already reach a distance of three times the half-mass radius from the cluster center before 1.5 Myr of evolution. Also the ejection rate is higher during the earlier times of evolution. Note that the time at which a star is classified as being ejected is different to when the ejection has actually happened. The time presented in the figure comprises the traveling time the stars need to reach to the ejection criterion and thus the ejections actually happened at an earlier time than shown in the figure. This can naturally explain the steeper increase of the number of ejected stars at a later time for the Banerjee et al. (2012) model in which the initial half-mass radius of the clusters is larger (${{r}_{h}}=0.8$ pc) than in the above two models (${{r}_{h}}=0.3$ pc). Thus an ejected star requires more time to reach the ejection criterion in their model (see Section 3 for the criteria). A rough estimate of the time when the ejections occur peaks at around 1 Myr, though it can vary with different initial conditions (S. Oh et al. 2015, in preparation).

Zoom In Zoom Out Reset image size

5.2. The Binary Fraction Among the Ejected O Star Systems

Not only single O stars but binaries containing O star component/s are ejected from star clusters via dynamical interactions. The binary fraction of ejected O stars can be substantial, up to $\approx 78$% (Table A1). Even initially single-star clusters eject O star binaries dynamically (Table A1). In Figure 11 is plotted the average binary fraction among ejected O star systems from the binary-rich cluster models as a function of cluster mass. A smaller fraction of binaries is ejected from the more massive clusters because the kick required to remove a system from the cluster is larger and more binaries are disrupted through interactions with other stars in the more massive clusters.

Zoom In Zoom Out Reset image size

The binary fraction among ejected O stars can be significantly altered by assuming a different initial period distribution. An initially higher fraction of close binaries leads to a higher binary fraction of ejected O stars. Thus the more realistic MS3OP_SP and MS3UQ_SP models produce higher binary fractions of ejected O star systems compared to the other models (Figure 11). The observed spectroscopic binary fraction among the Galactic field O stars is 43–50% (Gies 1987; Mason et al. 1998, 2009; Chini et al. 2012). This may be reproduced readily through dynamical ejections of O stars with a more realistic (than the Kroupa 1995b) period distribution function for O star binaries, such as constrained recently by Sana et al. (2012) (Equation (3)). This comes about because, as can be seen for the most realistic sequence (MS3UQ_SP), the dominant contributors to the population of ejected O star systems are $\leqslant {{10}^{4}}\;{{M}_{\odot }}$ clusters (lower panel of Figure 7), and these ejected O star systems have a binary fraction $\gtrsim 30$% (upper panel of Figure 11). It should be noted that the given observed binary fractions are lower limits of the binary fraction of the field O stars since only spectroscopic binaries are counted.

In the case of runaways, the binary fractions from our most realistic sequences, MS3OP_SP and MS3UQ_SP, are well consistent with the observed spectroscopic binary fraction of the runaway O stars in Gies (1987) and Mason et al. (1998, 2009), 19–29% (gray area in the lower panel of Figure 11), but too small compared to the value in Chini et al. (2012),³ 69% (gray horizontal line in the same figure). This may be due to the difference in assigning runaways between the Chini et al. (2012) work and the other studies. In Chini et al. (2012), stars are classified as runaways when they can be traced back to any cluster or associations in the study of Schilbach & Röser (2008). In other words, the runaways in their sample may be simply ejected O stars from star clusters which have masses $\lesssim {{10}^{4}}\;{{M}_{\odot }}$ since very massive clusters are exceedingly rare. The other studies assign stars to runaways with a large distance from the Galactic plane or with a high peculiar (radial) velocity.

In most of our cluster models (with the exception of the MS3OP_SP and MS3UQ_SP models) the binary fraction of O stars, whether they are inside a cluster or not, at 3 Myr is smaller than the observed O star binary fraction ($\approx 70$%, Chini et al. 2012; Sana et al. 2012, and references therein). This is the case because dynamical interactions between binaries and other cluster members disrupt binaries, preferentially those with a wide orbit (Kroupa 1995a; Parker et al. 2009; Marks et al. 2011) and because the Kroupa period distribution (which was constrained originally only for late-type stars) extrapolated in some of the model sequences to the O star regime has a larger fraction of long period binaries compared to the observed distribution of O star binaries (Figure 1).

While the mass of clusters which eject O stars most efficiently is the same for all models (${{M}_{{\rm ecl}}}\approx {{10}^{3.5}}\;{{M}_{\odot }}$), it has been shown that the dynamical ejection of O stars can be highly sensitive to the initial conditions of star clusters. In this section, we further discuss the effects of our choice of initial parameters for massive stars and for the clusters, and the available observational constraints.

The ejection fraction is strongly dependent on the initial cluster radius. As seen in Figure 2 the shapes of the curves look similar but the values of the ejection fraction are larger for initially smaller sized clusters. Thus it is important to check if our models are comparable to observed young star clusters. The half-mass densities of our N-body models at 3 Myr range from 175 to $\approx {{10}^{5}}\;{{M}_{\odot }}\;{\rm p}{{{\rm c}}^{-3}}$, from $\approx 160$ to $\approx 3.4\times {{10}^{4}}\;{{M}_{\odot }}\;{\rm p}{{{\rm c}}^{-3}}$, and from $\approx 50$ to $\approx 4000\;\;{{M}_{\odot }}\;{\rm p}{{{\rm c}}^{-3}}$ for initial ${{r}_{h}}=0.1$, 0.3, and 0.8 pc clusters, respectively. These values are well within the densities of observed young ($\lt 5$ Myr) massive (${{M}_{{\rm ecl}}}\gt 6000\;{{M}_{\odot }}$) star clusters, classified as star burst clusters in Pfalzner (2009), which range from 1500 to $4\times {{10}^{5}}\;{{M}_{\odot }}\;{\rm p}{{{\rm c}}^{-3}}$ (Figer 2008). However, caution is due in directly comparing the model cluster densities to the observed values since the observed densities are closer to the central (i.e., highest) densities rather than to the half-mass density. Concerning the initial radii of the embedded clusters, these are dynamically sub-pc (Testi et al. 1999; Lada & Lada 2003; Kroupa 2005; Smith et al. 2005; Tapia et al. 2009, 2011, 2014; Marks & Kroupa 2011).

The initial binary population plays a major role in the dynamical ejections of massive stars (see the differences between the MS3OP and MS3S models). Thus it is important to adopt realistic initial binary parameters to understand massive star ejections in reality. The pairing method for massive binaries adopted in most of our models results in extreme mass ratios, e.g., most of them are larger than 0.9. This was motivated by the high fraction of twins in observed massive binaries (Pinsonneault & Stanek 2006; Kobulnicky & Fryer 2007, but see also Lucy 2006), while recent observational studies favor a rather flat mass-ratio distribution for O-type binaries (Kiminki & Kobulnicky 2012; Sana et al. 2012) still excluding random pairing from the IMF for this spectral type. Adopting a more realistic initial mass-ratio distribution for O-type binaries may reduce the efficiency of the ejection, though the shape of the ${{f}_{{\rm ej},{\rm O}}}({{M}_{{\rm ecl}}})$ ejection fraction as a function of parent birth-cluster mass is expected to be kept. In this contribution we include a realistic model for the observed period distribution of O star binaries (Equation (3)). The model with the observed period distribution, which reproduces the observed abundance of short period binaries, exhibits a substantial increase of the ejection fraction, by $\approx 30$%, from massive clusters (${{M}_{{\rm ecl}}}\geqslant {{10}^{4}}\;{{M}_{\odot }}$). Whether the observed distribution is a real initial distribution is still unclear as the observed distribution has been compiled from several clusters which may have already evolved dynamically despite their young ages (1–4 Myr). The authors suggested that dynamical evolution would be negligible for such massive close binaries (${{{\rm log} }_{10}}(P/{\rm days})\lt 3.5$). However, the population of massive binaries can evolve through the dynamical interactions with other massive stars/binaries in the clusters (S. Oh et al. 2015, in preparation). Moreover, short-period massive binaries are vulnerable for merging either due to stellar evolution or the perturbation-induced collision of the two binary-star companions. Thus the effect of dynamical evolution needs to be taken into account iteratively to find the true initial period distribution of massive binaries. But this is beyond the aim of this study. For late-type binaries this has already been performed successfully by Kroupa (1995a, 1995b, see also Marks et al. 2011; Marks & Kroupa 2012).

Unlike the other initial conditions mentioned above, initial mass segregation does not result in a difference of the O star ejection fraction for clusters with ${{M}_{{\rm ecl}}}\geqslant {{10}^{3}}\;{{M}_{\odot }}$ and initially ${{r}_{h}}=0.3$ pc due to efficient dynamical mass segregation. But differences do occur for lowest mass clusters: initially mass-segregated clusters more efficiently eject O stars than non-segregated clusters, as dynamical mass segregation is inefficient to gather massive stars in the core of such clusters. For clusters with larger initial ${{r}_{h}}$ values, as the dynamical mass segregation time-scale becomes longer, the initial mass segregation can play a significant role. In observations, post gas-expulsion young massive clusters (age $\lesssim 5$ Myr) are generally compact (effective radii ≲ a few pc, Tapia et al. 2003; Kroupa 2005; Portegies Zwart et al. 2010).

Throughout this study the ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation is adopted. It follows from an intuitive theoretical concept that the mass of a star-forming molecular cloud core is finite and that star formation is a self-regulated process (Adams & Fatuzzo 1996), and it follows from the data on young stellar populations (Weidner et al. 2010, 2013a; Kirk & Myers 2011). The existence of a physical ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation is thus taken to be established (Weidner et al. 2014). Although the theory has been successful in describing the stellar contents of whole galaxies (e.g., Weidner et al. 2013b), it has been debated if the relation exists (see Section 2.1 in Weidner et al. 2013a for details) or if star formation is a purely stochastic process (e.g., Krumholz 2015). This is a fundamental issue for understanding galaxy-wide IMFs and how massive stars form, and we here thus stress that the data do not support the hypothesis that star formation is stochastic (see Section 7 of Kroupa 2015). For example, Hsu et al. (2012, 2013) presented that L1641, the low-density star-forming region (containing as many as $\approx 1600$ stars) of the Orion A cloud, is deficient of massive (O and early B) stars compared to the canonical IMF with 3–4σ significance and that with a probability of only 3% the southern region of L1461 and the ONC can be drawn from the same population, supporting that the high-mass end of the IMF is dependent on environmental density. In Section 4.2 of Krumholz (2015) three references (Calzetti et al. 2010; Fumagalli et al. 2011; Andrews et al. 2013) are cited as containing evidence against a physical ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation. Weidner et al. (2014) critically discuss these works, showing Andrews et al. (2013) data to be entirely consistent with the physical ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation. As already discussed in Section 3, whether the ${{m}_{{\rm max} }}$–${{M}_{{\rm ecl}}}$ relation is adopted or not does not affect the results obtained here though.

We do not include a gas potential, which would subsequently be removed rapidly, in the calculation (cf., Kroupa et al. 2001; Banerjee & Kroupa 2014). The main outcomes of gas expulsion are significant stellar loss from and expansion of clusters and their possible subsequent dissolution. Therefore gas expulsion can weaken the ejection efficiency by lowering the (central) density of a cluster if that occurs at a similar time when a majority of the ejection processes take place. In the future, further studies including a gas potential are needed to understand the effect of gas expulsion on the ejection of massive stars. Inclusion of a gas potential will be possible once the statistical quantification of the gas expulsion process is in place and computable (cf., Banerjee & Kroupa 2014). It has been argued recently that gas expulsion, if the process occurs at all, may have no effect on star cluster evolution (e.g., see Longmore et al. 2014, and references therein). It may, however, be difficult to detect the expanding signature resulting from residual gas expulsion since revirialization after gas expulsion can be rapid for massive clusters, ${{M}_{{\rm ecl}}}\gt {{10}^{3}}\;{{M}_{\odot }}$ (Portegies Zwart et al. 2010; Banerjee & Kroupa 2013).

To constrain which initial conditions are responsible for the observed field O star populations, one needs the observed quantities to compare with the results from models. However, it is very difficult to find which clusters individual field O stars originate from. Some ejected O stars can travel further than a few hundred pc away from their birth cluster within a few Myr, thus, without knowledge of their full-3D kinematics and the Galactic potential, it is almost impossible to find their birth cluster. Gaia - may help to find the birth clusters of the field O stars. However, the two-step ejection process of Pflamm-Altenburg & Kroupa (2010) complicates this. Their calculation suggests that 1–4.6% of O stars would appear to have formed in isolation because of the two-step ejection process. They used the observed values of the O star runaway fraction (10–25%, Gies & Bolton 1986; 46%, Stone 1991), and the binary fraction among the runaways (10%, Gies & Bolton 1986). Since most, if not all, ejected O star–O star binaries may go through the two-step ejection process and then cannot be traced back to their birth cluster, the fraction of O stars which appear to have formed in isolation can be higher if a binary fraction of the dynamically ejected stars is substantial. A quantitative study is required to determine the fraction of O stars experiencing the two-step ejection process.

We study how the ejection of O stars varies with cluster mass (or density as we assume a constant size for different cluster masses) under diverse initial conditions. We find that moderately massive clusters (${{10}^{3.5}}\;{{M}_{\odot }}$), which have formed about 10 O star systems, most efficiently shoot out their O stars to the field through energetic encounters between massive stars in the cluster core compared to other ones with other masses, for all model sequences, i.e., this result obtains independently of radius or the presence of initial mass segregation or binaries. Up to on average $\approx 25$% of the initial O star content is ejected from the clusters with this cluster mass in ${{r}_{h}}=0.3$ pc models (52% in ${{r}_{h}}=0.1$ pc clusters, Figure 2). However, the spread is large and ejections of all O stars leading to remnant young clusters void of O stars occur in about 1% of all such clusters (Figure 2). Fujii & Portegies Zwart (2011) suggested that one bully binary, dynamically formed, is mainly responsible for flinging out the stars from initially single-star clusters. Here we show that the number of ejections is significantly higher when the highly significant initial massive binary population is accounted for. The ejections, thus, occur via close encounters involving several binaries in reality since the massive stars are observed to have high multiplicity fractions in young star clusters (Chini et al. 2012; Sana et al. 2012, and references therein). Furthermore, more realistic clusters, which have O star binaries with periods following the Sana et al. (2012) distribution (MS3OP_SP and MS3UQ_SP), lead to even larger ejection fractions, especially for massive clusters. The mass-position of the peak of ${{f}_{{\rm ej},{\rm O}}}={{f}_{{\rm ej},{\rm O}}}({{M}_{{\rm ecl}}})$, however, does not change (${{M}_{{\rm ecl},{\rm peak}}}\approx {{10}^{3.5}}\;{{M}_{\odot }}$).

The decrease of the O star ejection fraction with increasing cluster mass for massive (${{M}_{{\rm ecl}}}\geqslant {{10}^{3.5}}\;{{M}_{\odot }}$) clusters can be represented by a simple power law of cluster mass with an exponent between −0.16 and −0.43. Clusters with an initial half-mass radius ${{r}_{h}}=0.3$ pc and with relatively simple initial conditions (Kroupa initial period distribution for O stars or no initial binaries) have an exponent of $\approx -0.4$, being in good agreement with −0.5 which is derived from the number of O stars from the IMF and the binary–single star scattering rate (Equation (12)). Models with a high fraction of close initial binaries or smaller ${{r}_{h}}$ deviate from the expected theoretical value by having larger ejection fractions.

We show that not only the ejection fraction but also the properties of the ejected O star systems depend on the initial cluster mass. The velocity distribution of the ejected O star systems shifts to higher velocities for more massive clusters, and so does their distance distribution from the cluster. But the distribution of ejected system masses varies weakly with cluster mass. The binary fraction of the ejected systems decreases with increasing cluster mass as does the binary fraction of the systems remaining in the cluster because binary disruption via interactions with other members of the cluster is more efficient for a denser cluster, i.e., more massive cluster. The binary fraction among the ejected O star systems is substantial, especially in the more realistic models in which the period distribution of the observed O star binaries is implemented (MS3OP_SP and MS3UQ_SP).

Considering that the cluster mass function is approximately a power law with a near Salpeter index, the main fraction of field O stars dynamically ejected from clusters originate from low/intermediate mass clusters (${{M}_{{\rm ecl}}}\lesssim {\rm a}$ few $1000\;{{M}_{\odot }}$, see Section 5). The dynamical ejection process populates $\approx 16$% to 38% of all O stars, formed in clusters, to the field. Combining our results on the dynamical ejection fraction and the binary-supernova scenario computed by Eldridge et al. (2011) it follows that the observed fractions of the Galactic field and runaway O stars are well accounted for theoretically. It should be noted that very massive ($\gtrsim 100\;{{M}_{\odot }}$) runaways can only originate from massive star-burst clusters (${{M}_{{\rm ecl}}}\gtrsim {{10}^{4}}\;{{M}_{\odot }}$) because only massive clusters can harbour a sufficient number of very massive stars to eject some of them.

We are grateful to Sverre Aarseth for making his NBODY6 code freely available and for continuing its improvements. We also thank Rolf Chini for comments on the observational data in his paper.

Here we append the table containing the results and properties of each model cluster (Table A1).

Table A1.  Results of all Models at 3 Myr

${{M}_{{\rm ecl}}}$ (${{M}_{\odot }}$)	${{m}_{{\rm max} }}$ (${{M}_{\odot }}$)	$\langle {{r}_{h}}\rangle$ (pc)	${{\langle }^{{\rm com}}}{{N}_{{\rm O}}}\rangle$	${{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle$	${{\langle }^{{\rm com}}}{{f}_{{\rm ej},{\rm O}}}\rangle$	$\langle {{f}_{b,{\rm O}}}\rangle$	$\langle {{f}_{b,{\rm ej},{\rm O}}}\rangle$	${{N}_{{\rm run}}}$
MS3OP
${{10}^{2.5}}$	21.2	0.60	1.11	0.09	0.080	0.84	0.78	100
${{10}^{3}}$	43.9	0.72	2.96	0.61	0.172	0.66	0.33	100
${{10}^{3.5}}$	79.2	0.70	9.57	2.38	0.252	0.45	0.17	100
${{10}^{4}}$	114.7	0.50	30.10	4.40	0.149	0.23	0.13	10
${{10}^{4.5}}$	136.2	0.50	105.00	11.25	0.105	0.15	0.05	4
NMS3OP
${{10}^{2.5}}$	21.2	0.56	1.14	0.03	0.020	0.87	0.50	100
${{10}^{3}}$	43.9	0.68	3.11	0.54	0.166	0.70	0.25	100
${{10}^{3.5}}$	79.2	0.65	9.26	2.21	0.244	0.51	0.15	100
${{10}^{4}}$	114.7	0.57	30.00	5.30	0.177	0.32	0.17	10
${{10}^{4.5}}$	136.2	0.50	98.25	9.25	0.092	0.27	0.12	4
MS3OP_SP
${{10}^{2.5}}$	21.2	0.62	1.23	0.09	0.090	0.79	0.67	100
${{10}^{3}}$	43.9	0.74	2.92	0.71	0.228	0.65	0.43	100
${{10}^{3.5}}$	79.2	0.72	9.11	2.32	0.257	0.51	0.32	100
${{10}^{4}}$	114.7	0.66	26.40	5.70	0.218	0.40	0.27	10
${{10}^{4.5}}$	136.2	0.55	96.50	14.00	0.145	0.27	0.21	4
MS3UQ_SP
${{10}^{2.5}}$	21.2	0.58	1.34	0.12	0.097	0.86	0.67	100
${{10}^{3}}$	43.9	0.73	3.32	0.62	0.177	0.66	0.56	100
${{10}^{3.5}}$	79.2	0.69	10.28	2.23	0.209	0.59	0.32	100
${{10}^{4}}$	114.7	0.63	33.60	6.30	0.188	0.45	0.29	10
${{10}^{4.5}}$	136.2	0.54	103.00	10.50	0.104	0.36	0.26	4
MS3S
${{10}^{2.5}}$	21.2	0.44	1.17	0.00	0.000	0.84	0.00	100
${{10}^{3}}$	43.9	0.55	2.88	0.17	0.043	0.56	0.13	100
${{10}^{3.5}}$	79.2	0.59	10.40	1.52	0.145	0.22	0.04	100
${{10}^{4}}$	114.7	0.52	34.30	3.30	0.096	0.06	0.00	10
${{10}^{4.5}}$	136.2	0.48	114.50	6.25	0.054	0.02	0.04	4
MS8OP
${{10}^{2.5}}$	21.2	0.93	1.10	0.00	0.000	0.97	0.00	100
${{10}^{3}}$	43.9	1.02	2.45	0.12	0.036	0.92	0.26	100
${{10}^{3.5}}$	79.2	1.05	8.81	0.49	0.053	0.54	0.20	100
${{10}^{4}}$	114.7	1.02	27.90	1.10	0.038	0.39	0.00	10
${{10}^{4.5}}$	136.2	1.00	86.50	2.25	0.026	0.35	0.00	4
${{10}^{5}}$ a	145.3	0.90	207.25	17.50	0.084	0.62	0.05	4
MS1OP
${{10}^{2.5}}$	21.2	0.60	1.16	0.28	0.248	0.91	0.64	100
${{10}^{3}}$	43.9	0.63	3.09	1.48	0.415	0.72	0.35	100
${{10}^{3.5}}$	79.2	0.49	9.10	4.55	0.516	0.47	0.15	100
${{10}^{4}}$	114.7	0.40	30.90	14.30	0.468	0.22	0.05	10
${{10}^{4.5}}$	136.2	0.33	106.00	33.75	0.320	0.12	0.04	4

Notes. Columns 1 and 2 show the initial stellar mass of the model clusters and the initial maximum stellar mass, respectively. Columns 3 and 4 present the averaged half-mass radius, $\langle {{r}_{h}}\rangle$, and the average number of O star systems present in the calculations, ${{\langle }^{{\rm com}}}{{N}_{{\rm O}}}\rangle$, at 3 Myr, respectively. Columns 5 and 6 are the average number of ejected O star systems, ${{\langle }^{{\rm com}}}{{N}_{{\rm ej},{\rm O}}}\rangle$, and the average O star system ejection fraction, ${{\langle }^{{\rm com}}}{{f}_{{\rm ej},{\rm O}}}\rangle$. The average binary fraction of O star systems remaining in the cluster, $\langle {{f}_{b,{\rm O}}}\rangle$, and the average binary fraction of ejected O star systems, $\langle {{f}_{b,{\rm ej},{\rm O}}}\rangle$, at 3 Myr are listed in columns 7 and 8, respectively. Note that the fraction of binaries among all O star systems is smaller than the initial value of ${{f}_{b,{\rm O}}}=1$ because systems are disrupted in these models. The number of realizations, ${{N}_{{\rm run}}}$, for each cluster mass is listed in the last column. Only clusters with ${{M}_{{\rm ecl}}}\geqslant {{10}^{2.5}}\;{{M}_{\odot }}$ are listed as less massive clusters do not have O stars in our models.

^aThe ${{10}^{5}}\;{{M}_{\odot }}$ cluster model is adopted from Banerjee et al. (2012). The initial conditions of the model from Banerjee et al. (2012) are slightly different from ours. Their O star binaries are mostly close binaries (see Figure 1) and this results in higher ejection fractions and fewer dynamical disruptions of O star binaries (higher binary fraction of O star systems at 3 Myr) compared to our models.

Download table as:  ASCIITypeset image