Runaway OB Stars in the Small Magellanic Cloud: Dynamical versus Supernova Ejections

Regarding the DES, there are many numerical studies in the literature. The main driver for generating massive runaways is the formation of a very dense central core of stars due to gravothermal collapse, which drives more massive stars toward the center (Fujii & Portegies Zwart 2011; Oh & Kroupa 2016). This greatly increases the probability of close encounters that can slingshot stars into the field. The cluster core collapse is halted only by the energy from hard binaries, either newly formed or primordial, which act as a kinetic energy source for the core (Fujii & Portegies Zwart 2014).

Fujii & Portegies Zwart (2011) conducted N-body simulations with a range of cluster masses from 6300 to 200,000 M_⊙, for a Salpeter initial mass function (IMF) and fixed central density. Their simulations do not include primordial binaries, although they do consider binaries that form dynamically. They find that stellar ejections are strongly dominated by a "bully" binary located in the core, that is relatively wide and composed of the most massive stars in the cluster, yielding the greatest interaction cross section. The bully binary is naturally produced via the gravothermal collapse of the cluster core, which occurs within about 1 Myr for the clusters simulated. Since the bully binary dominates the cross section for interactions, the encounter probability depends only weakly on the masses of the other, single stars. Therefore, the mass distribution of the runaway OB stars is not strongly modified; they find that the mass distribution for low-mass runaway stars <8 M_⊙ is consistent with the Salpeter slope. However, massive runaway stars >8 M_⊙ are significantly overrepresented, in particular for runaways >30 M_⊙. In addition, Fujii & Portegies Zwart (2011) examine the mean escape velocity as a function of mass for the runaways produced after 3 Myr, which gives more than enough time for the bully binary to form. They find a relatively constant median velocity near 42 km s⁻¹ for ejected stars >20 M_⊙.

Perets & Šubr (2012) carry out N-body simulations in which all massive stars reside in primordial binaries. They model a 5000 M_⊙ cluster with stellar masses ranging from 0.2 to 80 M_⊙ following a Salpeter IMF, and the binaries have initially zero eccentricity and semimajor axes set between 0.05 and 50 au. The cluster is evolved to an age of 2.7 Myr. They find that the velocity distribution of escaping stars having velocities in the range of 20−150 km s⁻¹ is independent of binary separations and cannot be produced via single−single encounters alone, further suggesting the presence of the bully binary that dominates the interaction cross section, in agreement with Fujii & Portegies Zwart (2011). Perets & Šubr (2012) find that the runaway fraction increases with mass, with the O star runaway fraction being two to three times greater than that of B stars, and the more massive stars also have higher ejection velocities. Furthermore, the DES can also generate runaway, noncompact OB binaries. The binary fraction of DES runaways decreases with ejection velocity, falling from ∼40% at velocities of 30 km s⁻¹ to ∼10% at 150 km s⁻¹.

Oh & Kroupa (2016) carry out comprehensive simulations of 10^3.5 M_⊙ clusters evolved to 3.5 Myr age, studying the effects of initial conditions, including mass segregation, binary fraction, period distribution, binary mass ratios, and eccentricities. They find that high stellar density is the dominant parameter driving high runaway frequencies, which is also aided by high primordial binary fractions, since ejections happen early, peaking around ages of ∼1 Myr. As also found by others (e.g., Fujii & Portegies Zwart 2011; Perets & Šubr 2012), O stars typically are ejected more frequently than B stars, therefore causing the mass function of all ejected OB stars, including runaways and walkaways, to be flatter than the parent IMF. However, the frequencies of only runaway O versus B stars may be more similar, depending on the mass ratios and period distribution; in particular, tight pairs of O star "twins" are the most dynamically active and produce the fastest runaways. In general, O stars tend to have faster runaway velocities, in some cases exceeding 200 km s⁻¹ (e.g., Perets & Šubr 2012). Oh & Kroupa (2016) note that the peak of the velocity distribution is closely related to the cluster mass and density, since the cluster potential determines the escape velocities. Thus, accounting for the cluster mass distribution would act to weight the velocity distribution toward lower values (Oh et al. 2015). The ejected binary frequencies are typically around ∼30% and are biased toward short periods.

The BSS is also well studied. For massive binaries, the core-collapse SN of the more evolved star results in one of two outcomes: (1) an unbound OB companion with an ejection velocity similar to its final pre-SN orbital velocity plus a poorly constrained SN kick, or (2) a bound binary consisting of the newly formed compact object and an OB companion, possibly observable as an HMXB. The reversal of the mass ratio prior to the explosion and widening of the orbit tend to inhibit acceleration to runaway velocities (v > 30 km s⁻¹). Brandt & Podsiadlowski (1995) carried out simulations of massive binary systems to evaluate the post-SN binary and HMXB properties. They modeled the outcomes for a binary system of fixed initial stellar masses and a constant logarithmic distribution of initial orbital periods, with an SN kick velocity distribution based on the observed pulsar birth velocities. They demonstrate that the period anticorrelates with both kick velocity and kick directions opposing the orbital motion, resulting in a strong anticorrelation between the final orbital separation and post-SN systemic velocity. The bulk of their massive binaries that remain bound have periods of <100 days, which also is roughly the threshold for 30 km s⁻¹ runaway velocities. About 1/4 of these systems obtain eccentricities leading to strong dynamical mass transfer, inducing merging or disruption; thus, ultimately about 20% of their massive binaries remain stably bound. They also note that the remaining star's spin constitutes a substantial fraction of the total orbital angular momentum, implying high rotation velocities.

More recently, Renzo et al. (2019) perform comprehensive binary population synthesis modeling to study populations of BSS products, including unbound single stars. The binaries are mainly characterized by three parameters: (1) the primary-star masses, ranging from 7.5 to 100 M_⊙ and weighted by a Kroupa IMF; (2) a flat distribution in initial mass ratio q, ranging from 0.1 to 1; and (3) an empirically motivated, −0.55 power-law distribution for the initial orbital period over a range of log(0.15/days) to log(5.5/days). In the fiducial model they assume a metallicity of Z = 0.02; however, they also consider lower metallicities. They also treat a variety of evolutionary processes, including wind mass loss, mass-transfer efficiency, common envelope evolution, and SN kick parameters.

The results of Renzo et al. (2019) are generally consistent with those of Brandt & Podsiadlowski (1995) for bound runaways. For the unbound OB stars, they find that in general the ejection velocities are too slow to produce large numbers of runaways with velocities >30 km s⁻¹. The mean ejection velocity for their fiducial population of OB companions is ∼12 km s⁻¹, and toward higher companion masses, the mean ejection velocities are even lower. This is indeed the opposite trend that we expect from the DES as discussed above. This robust prediction for low runaway frequencies is partly due to the fact that essentially all massive runaways have experienced mass transfer such that the secondary becomes more massive than the primary, widening the orbit before the primary explodes. Thus, both the increased mass of the ejected star and the lower orbital velocity contribute to the slower ejection speeds. In their fiducial model at solar metallicity, 67% of the OB companions are disrupted, such that only 2.5% of the parent binaries are ultimately ejected as single runaways with v > 30 km s⁻¹ and 48% have v < 30 km s⁻¹. These simulations give an 11% frequency for post-SN bound systems, out of all OB binaries. The remaining 39% of the parent binaries result in stellar mergers pre-SN (22%) and evolved single runaways and walkaways (17%). Kochanek et al. (2019) also carried out binary population synthesis models that are largely consistent with these branching ratios.

We will use these predicted kinematic trends for the DES and BSS to discriminate between these two ejection mechanisms. A third, hybrid "two-step" mechanism also generates runaways that correspond to a small subset of both DES and BSS populations (Pflamm-Altenburg & Kroupa 2010). We will also evaluate the effect of these objects on the relative total contributions of each mechanism.

We expect the highest-velocity stars to be due to dynamical ejections (Section 1.1). Based on their models described above, Perets & Šubr (2012) predict a runaway velocity distribution of v^−3/2 below ∼150 km s⁻¹, which projected to 2D corresponds to ${v}_{\mathrm{loc}}$ of ∼122 km s⁻¹, steepening to v^−8/3 at higher velocities. Oh & Kroupa (2016) find similar slopes of −1.4 to −2.1, fitted over the entire runaway range. These works do not consider the effects of steady-state star formation with a cluster mass function, but Oh & Kroupa (2016) find that O- and B-star velocity distributions are similar to each other at the highest velocities, and so the composite velocity distribution for runaways should be quite robust. In Figure 1, we see that the predicted velocity relation for ${v}_{\mathrm{loc}}$ < 122 km s⁻¹ from Perets & Šubr (2012) agrees well with our runaway velocity distribution, although we note that the observed distribution is somewhat steeper, as would be expected by significant contamination from BSS runaways at lower velocities (see Sections 2.2 and 4.2). Thus, the overall distribution is consistent with our finding in Paper I that the DES mechanism dominates the runaway population in our sample.

Zoom In Zoom Out Reset image size

Adopting the predicted ${v}_{\mathrm{loc}}$^−3/2 distribution below 122 km s⁻¹, we can therefore roughly estimate the total number of dynamically ejected runaways out of our sample of 304 field OB stars by extrapolating from the high-velocity tail. Assuming that the 40 stars with 75 km s⁻¹ ≤ ${v}_{\mathrm{loc}}$ < 122 km s⁻¹ are all dynamically ejected runaways, we obtain 108 DES runaways in the range 30 km s⁻¹ ≤ ${v}_{\mathrm{loc}}$ < 75 km s⁻¹. There are a total of 149 stars observed in this velocity range, therefore implying that 41 stars, or 21% ± 7%, are due to the BSS. Adding the 40 stars with 75 km s⁻¹ ≤ ${v}_{\mathrm{loc}}$ < 122 km s⁻¹ and also the 8 highest-velocity stars (${v}_{\mathrm{loc}}$ > 122 km s⁻¹), which are all assumed to be dynamically ejected, to the 108 predicted DES runaways in the range 30 km s⁻¹ ≤ ${v}_{\mathrm{loc}}$ < 75 km s⁻¹, we obtain a total DES contribution of 79% ± 8% (=156/197) to the runaway population. Runaways compose 65% ± 6% (=197/304) of our sample, and so DES runaways alone, without walkaways, are about half (51% ± 7%) of all our SMC field OB stars. The errors quoted are for Poisson statistics only and do not account for systematic uncertainties; for example, it remains possible that a few objects in the 75–122 km s⁻¹ velocity range are BSS runaways, implying a possible slight overestimate in our total DES frequency. One HMXB has ${v}_{\mathrm{loc}}$ in this range, supporting this possibility (Section 5.2).

We can also make an independent estimate of the contribution of BSS runaways to our sample. In a steady-state population with a constant star formation rate, the OB runaway fraction from BSS is smaller than that from DES. Renzo et al. (2019) find that 1.2% and 10% of OB stars with m > 15 M_⊙ in a steady-state population at SMC metallicity are ejected single-star runaways and walkaways, respectively. However, the post-SN bound frequency is 13% for companions to black holes (BHs) or neutron stars (NSs), a large fraction of which achieve v > 30 km s⁻¹, depending on the assumed kick velocities and BH fallback prescriptions. Thus, a total of roughly 24% of OB stars in a steady-state population are the survivors of post-SN binary systems, and the remaining ∼76% are primaries that have not yet exploded. This is consistent with the findings of Moe & Di Stefano (2017), who estimated that ∼20% of OB stars in a steady-state population are the secondaries in which the true primaries have exploded as SNe.

About half of the post-SN systems (13%) remain bound. The models at SMC metallicity predict a total walkaway-to-runaway ratio of 11.1, so we therefore expect that ∼2% of all OB stars are BSS runaways. The RIOTS4 field stars correspond to about 25% of the SMC population (Oey et al. 2004), and therefore we expect that about 8% of our field stars are runaways due to the BSS. Similarly, since 65% ± 6% of our field stars are runaways, then BSS runaways are about 12% ± 2% of runaways. This is slightly lower than our estimate of 21% ± 7% based on the prediction for DES runaways (Section 2.1), and moreover, this breakdown does not account for objects that experience both DES and BSS runaway acceleration, including two-step runaways. We carry out a full accounting of the DES and BSS runaways below in Section 5.

Altogether, these first-order estimates suggest that the DES:BSS allocations are roughly in the range 80:20, respectively, among the SMC runaway OB stars. These values will be adjusted when we fully evaluate the ratio of DES to BSS runaways in Section 5.3.

It is possible for the DES to eject massive binaries, in addition to single stars. Therefore, noncompact binary runaways, such as EBs and SB2s, are a direct probe of the DES mechanism. Faster runaway binaries can be produced via binary−binary interactions in a dense cluster core, whereas binary interaction with a single star causes the binary typically to be ejected at much lower velocities (Fujii & Portegies Zwart 2011; Perets & Šubr 2012). We are thus in a position to compare the characteristics of our noncompact, massive binary runaways to predictions, to better understand their production and interaction histories.

Perets & Šubr (2012) find a maximum velocity of ≲200 km s⁻¹ for dynamically ejected binary runaways in their simulations, while some models of Oh & Kroupa (2016) have maximum speeds that are double this value. Taken at face value, our highest-velocity binary is an EB (M2002-81258) with ${v}_{\mathrm{loc}}$ = 201 ± 36 km s⁻¹ and a period of 2.7 days (Pawlak et al. 2016). Two more high-velocity binaries are an SB2 (M2002-36213) with ${v}_{\mathrm{loc}}$ = 121 ± 32 km s⁻¹ and an EB/SB2 star (M2002-65355) with ${v}_{\mathrm{loc}}$ = 109 ± 34 km s⁻¹ and period of 1.2 days (Pawlak et al. 2016). These high-velocity noncompact binaries are likely the result of a close interaction with a hard and/or massive bully binary in the cluster core. Although we caution that individual Gaia PMs may have unknown errors, the short periods are consistent with expectation for high-speed binaries. We find that, in general, the binary frequency of our EBs and SB2s decreases with increasing velocity, a result consistent with the simulations of Perets & Šubr (2012).

Among our 197 runaway stars with ${v}_{\mathrm{loc}}$ > 30 km s⁻¹, 3 are EBs, 5 are SB2s, and 3 are both EBs and SB2s. This yields a DES binary runaway frequency of >6% ± 2%, which is a substantial lower limit since there are likely many additional binaries we are unable to identify, not only because of selection effects but also because of the post-SN, two-step ejection mechanism (Section 5.2). Mason et al. (2009) give a runaway binary fraction of 43% and total field binary fraction of 59%, based on a comprehensive accounting of the observed field O stars. Lamb et al. (2016) also estimate a field OB binarity of 59% ± 12% for a small subset of the RIOTS4 survey, and Chini et al. (2012) find binary frequencies of 69% ± 11% and 43% ± 13% for O star runaways and field stars, respectively.

On the other hand, the predicted binary frequency of runaway O stars is ∼30% (Perets & Šubr 2012; Oh & Kroupa 2016), and slightly higher for lower-mass primaries. While this is consistent with our lower limit on the DES binary runaway frequency of 6%, the more comprehensive accounting in the above studies shows values that are roughly double the predicted ones. In any case, given that noncompact runaways must originate from dynamical ejections, the kinematics and frequency of this population are consistent with a dominant population of DES ejections, as found above.

Since the dynamical and SN ejection mechanisms predict contrasting relationships between mass and runaway velocity (Section 1.1), we obtain spectroscopically determined masses of the RIOTS4 stars to further evaluate the allocation between the two ejection processes. Effective temperatures (T_eff), luminosities (L), and stellar masses are calculated following an approach similar to that described in Lamb et al. (2013), but updated using the stellar evolutionary models for rotating stars of Brott et al. (2011). In particular, T_eff are derived according to the spectral types published in Lamb et al. (2016) with conversions of spectral type to T_eff from Massey et al. (2005) for O-type stars and from Crowther (1997) for B-type objects. Bolometric magnitudes (M_bol) and luminosities were estimated following Massey (2002), based on V magnitude and adopting a distance modulus (DM) of 18.9 (Harries et al. 2003). The bolometric correction (BC) was obtained as BC = 27.99–6.9 log T_eff (Massey et al. 2005). We use the extinction A_V listed in Lamb et al. (2013), extracted from the SMC extinction maps in the Magellanic Clouds Photometric Survey (Zaritsky et al. 2002). Stellar masses are estimated from the positions of the stars in the Hertzsprung–Russell (H-R) diagram, compared with the rotating (150 km s⁻¹), single-star, evolutionary tracks of Brott et al. (2011) for SMC metallicity. The masses are obtained by interpolating between the tracks for the evolved stellar mass values at the observed positions in the H-R diagram. For two stars, M2002-38024 and 59319, we adopt the masses from Lamb et al. (2013), since their positions relative to the Brott et al. (2011) evolutionary tracks do not permit reliable mass determinations.

Table 1 gives our mass estimates. The median ratio of our revised masses to those determined by Lamb et al. (2013) is 1.03 for the 107 stars in common, with standard deviation of 0.17. The greatest uncertainty in determining the stellar masses are the spectral classifications. Of our 297 stars for which we obtained masses, 238 have a well-determined spectral type and thus a well-constrained mass, while 59 do not. For these 59 stars, we adopt the average of the lower and upper mass limits obtained from the limits in the star's range of possible spectral type. The uncertain masses for these stars are flagged with a ":" in Table 1; these are mostly Be stars that lack a spectral type and SB2s, as well as a few Oe stars and "O + B" binaries. For the Be stars without a well-constrained spectral type, we calculated the average T_eff from adopting B0e and B2e spectral types. For SB2s, the masses are calculated for the hotter, more massive star; in the case of "twins," the luminosity is reduced by a factor of two and the quoted mass value is the average of these two. Masses for EBs not identified as SB2s are likely slightly overestimated since these often may be twins. For the four B[e] stars, we obtain spectral types from Graus et al. (2012). Errors for the 238 reliable spectral types are obtained by assuming half a spectral type on either side of that observed, and these errors may therefore be somewhat overestimated. Figure 2 displays the position of the stars in the H-R diagram, together with the Brott et al. (2011) evolutionary tracks.

Zoom In Zoom Out Reset image size

Table 1.  Kinematic Data and Fundamental Parameters for RIOTS4 Field OB Stars

IDa	SpTypeb	Subgroupc	${v}_{\mathrm{loc}}$ d	Error	Me	Errorf	${v}_{r}\sin \ i$	Error	Teff	Error	log L	Error
			(km s−1)	(km s−1)	(M⊙)	(M⊙)	(km s−1)	(km s−1)	(kK)	(kK)	log(L/L⊙)	log(L/L⊙)
107	Be3	−, −, −, e	20	21	14.9:	6.3	⋯	⋯	22.2	7.4	4.52	0.42
1037	B0.5 V	−, −, −, −	99	30	14.6	2.5	92	35	26.7	3.1	4.55	0.14
1600	O8.5 V	E, −, −, −	43	30	26.8	1.5	91	17	35.5	1.5	5.11	0.05
1631	B1e2	−, −, −, e	51	23	14.6	3.2	197	11	23.4	3.7	4.71	0.19
1830	B0.5 III	−, −, −, −	32	26	22.3	4.4	86	12	26.7	3.1	5.11	0.14
2034	B	−, −, −, −	21	26	17.2:	7.4	⋯	⋯	22.2	7.4	4.67	0.42
2093	B1e3+	−, −, −, e	63	27	13.1	2.8	⋯	⋯	23.4	3.7	4.52	0.19
3224	B1e2+	−, −, −, e	9	21	19.2	4.3	⋯	⋯	23.4	3.7	5.0	0.19
3459	O9.5 I	−, −, −, −	7	27	33.3	2.5	238	17	31.9	2.1	5.5	0.08
3815	Be2	−, −, −, e	147	29	17.2:	7.4	⋯	⋯	22.2	7.4	4.67	0.42

Notes.

^aFrom Massey (2002). ^bSpectral types are from Lamb et al. (2016), except for the four B[e] stars (M2002-29267, 46398, 62661, 83480), which are taken from Graus et al. (2012). ^c"E," "S," and "X" indicate EB, SB2, and HMXB, respectively; "e" indicates emission-line star (Oe or Be). ^dLocal residual transverse velocity from Paper I. ^eMasses that are uncertain are flagged with ":" (see Section 3.1). Masses flagged with "n" are one of the following: The mass for SMC X-1 (M2002-77458) is taken from Rawls et al. (2011), and those for M2002-38024 and 59319 are taken from Lamb et al. (2013). For the mass of the "mid"-Oe star M2002-73795 we calculated the average T_eff from adopting O5e and O7e spectral types, and for the Oe star M2002-75689 we adopted O3e and O9e spectral types. For Be stars without constrained spectral types, we calculate the average T_eff for B0e and B2e types. For the "O + B" binaries (M2002-11238, 22178, 66160) we calculated the average T_eff from adopting O3 and O9 spectral types, in estimating masses. ^fErrors are based on ± half a subtype. Errors for stars flagged with ":" are computed from the upper/lower limits on the spectral type range (see Section 3.1).

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

We caution that, in what follows, our analyses are based on masses derived from modeled stellar evolution, and not dynamical mass determinations. Moreover, stars in binary systems will have more uncertain masses, as noted in Table 1, and if they have experienced binary mass transfer, their positions on the H-R diagram may not provide a meaningful mass estimate (e.g., Wang et al. 2020). This is especially the case if an active accretion disk is present. In particular, we note that the well-known HMXB SMC X-1 has a well-determined dynamical mass of 15.35 ± 1.53 M_⊙ for the primary star, based on accurate X-ray eclipse data (Rawls et al. 2011), whereas our mass obtained from the H-R diagram is twice this value, 32.2 ± 4.5 M_⊙. This may be due to the accretion disk amplifying the observed T_eff and/or luminosity. Furthermore, the star is likely to be rejuvenated by binary mass transfer before the SN event (e.g., Vinciguerra et al. 2020), altering its stellar structure. Thus, we caution that the masses we obtained for the rest of our HMXBs might also be significantly overestimated.

A power-law fit to the stellar mass distribution for m > 20 M_⊙ gives a slope of −2.96 ± 0.34, which agrees with the present-day mass function (PDMF) slope of −2.8 to −3.1 for m > 20 M_⊙ obtained for the subset of RIOTS4 stars studied by Lamb et al. (2013).

Our stellar masses are plotted against ${v}_{\mathrm{loc}}$ in Figure 3, with the HMXBs, EBs, and SB2s identified. Four stars are both EBs and SB2s. As described in Section 1.1, the BSS runaway and walkaway velocities are expected to decrease with mass (Renzo et al. 2019), while for the DES the ejection velocity and runaway fraction of stars both increase with mass (Banerjee et al. 2012; Perets & Šubr 2012; Oh & Kroupa 2016). Therefore, we might expect (1) a slow-moving population, with a range of masses, but predominantly at low masses, resulting from the SN ejection scenario; and (2) a higher-mass, high-velocity population resulting from the DES. The density map in Figure 4 suggests the existence of two such populations, although they are clearly intermixed. We note that the slow population appears to be centered more densely near (30 km s⁻¹, 18 M_⊙); this corresponds roughly to our detection and completeness limits, while a second population is also suggested at higher mass and velocity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The mass functions for the nonrunaway population (${v}_{\mathrm{loc}}$ < 30 km s⁻¹) and the fast runaways (${v}_{\mathrm{loc}}$ > 75 km s⁻¹) are shown in Figure 5. Power-law fits to each stellar mass distribution for m > 20 M_⊙ are shown; the slope for the nonrunaway population is −2.95 ± 0.46, and the slope for the fast runaways is −2.13 ± 0.79. We note that the slope for the nonrunaway population is essentially identical to that of the total mass distribution for m > 20 M_⊙ (see above), while the slope for the fast runaways is considerably flatter. Although the error is large, this further supports our understanding that the fastest runaways are produced via the DES mechanism, which favors the ejection of more massive stars, resulting in a flatter mass distribution (Section 1.1).

Zoom In Zoom Out Reset image size

For the DES, Fujii & Portegies Zwart (2011) predict that the median velocity for ejected stars more massive than ∼20 M_⊙ should be relatively constant, a result also found by Banerjee et al. (2012), and one that contrasts with expectations for the BSS. Considering only our stars with well-constrained masses, for the mass ranges 18–25 M_⊙ (79 stars), 25–30 M_⊙ (47 stars), and >30 M_⊙ (41 stars), we obtain median velocities of 43, 40, and 42 km s⁻¹, respectively. Thus, the constant median velocity is again consistent with DES dominating our runaway population. We note that these median velocities are somewhat greater than those predicted by Fujii & Portegies Zwart (2011), which when projected to 2D are on the order of 34 km s⁻¹. This difference may be caused by our sample selection bias against walkaways. However, given the expected significant contamination of slower, BSS walkaways, it may also point to dominant origins in more massive clusters (e.g., Oh et al. 2015) or differences in, e.g., the primordial binary parameters assumed by Fujii & Portegies Zwart (2011) relative to those in the SMC.

It is interesting to compare our runaway population in Figure 3 with the models of Oh & Kroupa (2016). Their Figure A1 shows that the high-mass, high-velocity parameter space is sensitive to a variety of parameters. At face value, our data are more consistent with binary mass pairing corresponding to either a uniform mass ratio distribution or ordered pairing. Since dynamical processing in the cluster core likely disrupts most wide binaries, the DES preferentially produces runaway OB stars with close companions, which follow a uniform mass ratio distribution (Sana et al. 2012), whereas wide visual companions are weighted toward smaller mass ratios (Moe & Di Stefano 2017). This scenario provides enough massive binaries with fairly uniform mass ratios to generate the runaways with velocities on the order of 100–200 km s⁻¹. We also see noncompact binary runaways at both high mass and high velocity, which are also consistent with these parameters. However, we caution that the models are based on simulations of single-mass clusters, whereas the observed SMC runaways are products from many clusters of different masses, which will affect the velocity distribution (Oh et al. 2015).

Although, as discussed above, the runaways are dominated by the DES mechanism, we see in Figure 3 that 7 of the 14 HMXBs have ${v}_{\mathrm{loc}}$ > 30 km s⁻¹. HMXBs are direct tracers of the SN mechanism and, in particular, systems that remain bound, which have faster velocities than single BSS stars. Those with OB secondaries are predicted to have the highest BSS velocities, which can reach speeds of 80 km s⁻¹ or more (Brandt & Podsiadlowski 1995; Renzo et al. 2019). The median velocity ${v}_{\mathrm{loc}}$ of our 14 HMXBs is 28 km s⁻¹, which when projected to 3D is 34 km s⁻¹, consistent with Coe (2005), who finds an average space velocity of SMC HMXBs ≳ 30 km s⁻¹. This is greater than the predicted median systemic velocity for NS+MS binaries of only 20 km s⁻¹ found by Renzo et al. (2019). It is likely due in part to a selection effect favoring the fastest HMXBs, since the fastest runaways are the tightest (Brandt & Podsiadlowski 1995), and therefore most luminous, HMXBs. There are likely many undetected, bound BSS runaways and walkaways.

However, as we shall show in Sections 5.2 and 5.3, it is likely that the two-step ejection mechanism is responsible for a significant number of post-SN walkaways and runaways, including HMXBs. The expected ratio of BSS walkaways to runaways from the models of Renzo et al. (2019) is several times higher than our observed ratio of unity. Although our sample is biased against walkaways, the difficulty in generating runaway velocities by the BSS mechanism alone suggests that the two-step mechanism is important in producing our observed BSS runaway population.

The distribution of ${v}_{r}\sin \ i$ for OB stars is known to be bimodal, especially for B stars (e.g., Wolff et al. 2007; Bragança et al. 2012; Dufton et al. 2013), and we therefore expect to find such bimodality in our sample. We also expect higher average rotation speeds than typically found in the Milky Way, due to the low metallicity of the SMC, which suppresses a star's ability to lose angular momentum through stellar winds and therefore drives up rotational velocities (e.g., Maeder & Meynet 2000). The high expected rotation speeds and low metallicity also lead to a high frequency of classical Oe and Be stars, which spin much faster than non-Oe/Be stars (e.g., Rivinius et al. 2013; Golden-Marx et al. 2016). Oe/Be stars compose 40% ± 4% (=123/304 stars) of our field OB star sample.

The distribution of our measured ${v}_{r}\sin \ i$ is shown in Figure 6, for 140 non-Oe/Be stars (top) and 53 Oe/Be stars (bottom). O/Oe stars are shown in the left panels and B/Be stars in the right panels. Separating the 47 Be stars from the normal B stars, the bimodality in the ${v}_{r}\sin \ i$ distribution apparently corresponds to the different ${v}_{r}\sin \ i$ distributions for Be versus non-Be stars. Non-Be stars show a strong peak at low ${v}_{r}\sin \ i$, with decreasing numbers at higher values. In contrast, Be stars show a peak at much higher ${v}_{r}\sin \ i$, consistent with a population dominated by high spin velocities around 250–350 km s⁻¹ that is broadened by ${v}_{r}\sin \ i$ projection to lower apparent values. These are likely to be substantial underestimates due to gravity darkening of these highly oblate stars. The O stars are much fewer, but K-S tests show that the ${v}_{r}\sin \ i$ distributions agree statistically with those for the B stars. There are only six Oe stars with reliable ${v}_{r}\sin \ i$ measurements, due to their tendency to show He I emission, which often causes infill of the features used to measure ${v}_{r}\sin \ i$ (Golden-Marx et al. 2016).

Zoom In Zoom Out Reset image size

We can compare our ${v}_{r}\sin \ i$ distributions of isolated SMC OB stars to those of OB stars in the SMC Wing published by Ramachandran et al. (2019). They, too, find that Be stars are much faster rotators than normal OB stars. Our low-velocity peak is situated at our effective resolution limit of ∼80 km s⁻¹, lower than the low-velocity peak of ∼120 km s⁻¹ quoted by Ramachandran et al. (2019). They also report a mean ${v}_{r}\sin \ i$ of 230 km s⁻¹ for Be stars, which is higher than our value of 190 km s⁻¹. Our fastest-rotating star has ${v}_{r}\sin \ i$ of 439 km s⁻¹, substantially slower than their fastest star at ∼550 km s⁻¹. Thus, in general, their ${v}_{r}\sin \ i$ values are slightly higher than ours. Ramachandran et al. (2019) also use the iacob_broad code to determine ${v}_{r}\sin \ i$, but their sample includes stars of much later spectral types, extending to late B stars. The variation between our results suggests that lower-mass B stars may rotate somewhat faster than massive OB stars. This may be consistent with the later spectral types having weaker winds, causing them to retain angular momentum.

Stellar rotation provides another important parameter for discriminating between dynamical and SN ejections of OB runaways. As discussed above in Section 2, the vast majority of BSS ejections are expected to result in slower space velocities, and only 10%–20% of them are predicted to have runaway velocities ≳30 km s⁻¹ at SMC metallicity (Section 2.2). The runaway BSS ejections are those originating from the tightest interacting binaries, which have greater orbital velocities at the time of the SN explosion. This implies that mass transfer from the higher-mass primary will spin up the secondary star to rotation speeds near breakup velocity (e.g., Brandt & Podsiadlowski 1995; de Mink et al. 2009). Thus, after the primary explodes, the secondary is ejected with both high space and rotation velocities. This applies to both bound and unbound BSS runaways and is supported by observations (Blaauw 1993; Hoogerwerf et al. 2001; Walborn et al. 2014; Maíz Apellániz et al. 2018). In contrast, single-star dynamical ejections are less likely to be from tight binaries, and hence their rotation speeds should be similar to those of nonrunaway stars.

Figure 7 shows ${v}_{r}\sin \ i$ as a function of ${v}_{\mathrm{loc}}$ for the 201 RIOTS4 stars that have measurements of both quantities. There is a prominent population of stars with low ${v}_{r}\sin \ i$, near the IMACS resolution limit of 80 km s⁻¹, that have runaway velocities extending to high values. This population is strongly inconsistent with the SN ejection mechanism and must correspond to runaways accelerated by the DES. There are 52 stars with ${v}_{r}\sin \ i$ ≤ 100 km s⁻¹ out of the 128 runaway stars with ${v}_{\mathrm{loc}}$ ≥ 30 km s⁻¹, which sets a lower limit of 41% on the contribution of DES runaways, consistent with our inference that these strongly dominate over the contribution of BSS runaways (Section 2; Paper I).

Zoom In Zoom Out Reset image size

We can also use the Oe/Be stars to probe the role of SN ejections in generating runaways. Since Oe/Be stars are expected to be near breakup v_r, it seems likely that runaway Oe/Be stars have experienced mass transfer, and therefore many, if not most, acquired their high ${v}_{\mathrm{loc}}$ from SN kicks (e.g., de Mink et al. 2013; Boubert & Evans 2018). In Figure 7, Oe/Be stars are indicated with squares and non-Oe/Be stars with circles. Figure 7 shows that while most Oe/Be stars have ${v}_{\mathrm{loc}}$ < 75 km s⁻¹, many of them are in the runaway regime >30 km s⁻¹. On the other hand, all of the extreme runaways with ${v}_{\mathrm{loc}}$ > 90 km s⁻¹ are non-Oe/Be stars, supporting the expectation that only DES can produce the fastest runaways, and that these tend to be single. The connection between Oe/Be and the BSS is further supported by the fact that very few of our noncompact binaries, which are tracers of the DES mechanism, are Oe/Be stars: only 1 of our 15 EBs and 2 of our 11 SB2s are Oe/Be stars, while none of our runaway EBs and only one of our runaway SB2s are Oe/Be stars.

The tendency for fast rotators to be runaways is further shown by the fact that 33 out of our 57 stars with ${v}_{r}\sin \ i$ > 200 km s⁻¹ are runaways, or 58% ± 13%. Using their Galactic OB star sample, Maíz Apellániz et al. (2018) found that 13% ± 4% of their stars with ${v}_{r}\sin \ i$ > 200 km s⁻¹ are runaways, using a similar definition of runaways based on velocity dispersion. Since their sample includes both field and nonfield stars, the large difference between these two fractions tells us that OB stars with ${v}_{r}\sin \ i$ > 200 km s⁻¹ are much more likely to be runaways in the field than nonrunaways in clusters. This is consistent with our conclusion that most stars with large ${v}_{r}\sin \ i$ achieved such high rotation rates due to mass transfer in close binaries, which are more likely to produce runaways. However, we caution that only a subsample of our stars have measured ${v}_{r}\sin \ i$ available. Also, high ${v}_{r}\sin \ i$ is associated with low metallicity, as discussed in Section 4.1, and so a large fraction of DES runaways in our sample will also have high ${v}_{r}\sin \ i$.

Thus, we may crudely expect that most Oe/Be stars are post-SN runaways (e.g., McSwain & Gies 2005; de Mink et al. 2013; Boubert & Evans 2018). On the one hand, the frequency of Oe/Be stars underestimates that of the BSS runaways because not all SN ejections end up as Oe/Be stars, and moreover, the Oe/Be phenomenon likely has a finite lifetime, even for those BSS objects that become Oe/Be stars. On the other hand, runaway DES binaries are also most likely to be tight binaries that undergo mass transfer, possibly producing Oe/Be stars; and single Oe/Be stars may also spin up through other mechanisms, although these should be mostly late-type B stars, which can retain most of their angular momentum owing to weaker winds (e.g., Ekström et al. 2008; de Mink et al. 2009). The relative magnitude of these effects is not known, but they counteract each other. Thus, it may not be unreasonable to crudely assume that all runaway Oe/Be stars are accelerated by SNe, including two-step ejections, and the remaining runaways are due to dynamical ejections. As noted above, about half of Oe/Be stars have runaway velocities, as we also find for HMXBs (Section 3.2). This consistency further supports adopting Oe/Be stars as tracers of BSS ejections, and it again reveals an unexpectedly large number of BSS runaways.

There are 69 Oe/Be stars out of the 197 runaway stars in our full sample, therefore implying a BSS runaway fraction of 35% ± 5%. This crude estimate is larger than our ∼20% estimate obtained in Section 2. It does not account for objects in common between DES and BSS but suggests a slightly larger frequency (see Section 5.3).

In Section 2, we found that about 70% of our field OB stars are runaways, implying an overall frequency of 0.18 for a field OB star fraction of about 25% (Table 2, Total observed). We therefore take walkaways to be 30% of the field, but as noted in Section 2, we have a large selection bias against walkaways, requiring a correction factor of ∼2.4. This therefore implies that our observationally derived walkaway frequency is also 0.18, or roughly equal to that of the runaways (Table 2, Total observed). This uniform ratio of walkaways to runaways is harder to match with predictions, since both DES and BSS predict many more walkaways than runaways. One possible issue could be that the observed ratio is affected by spurious runaways owing to poorly quantified effects with Gaia PMs.

It is also plausible that there is further incompleteness among the walkaways, in addition to the selection bias. An additional observational bias against both walkaways and cluster stars is quite likely in the Oey et al. (2004) survey from which RIOTS4 is drawn, due to extinction and crowding near clusters (e.g., Schoettler et al. 2020). We note that the Oey et al. (2004) survey has a total of 1360 OB candidates, whereas, in contrast, the survey by Evans et al. (2004) identified about 2500 OB stars, including 23 O stars that were fainter than the RIOTS4 selection criterion. A large fraction of their OB stars have uncertain spectral types and are unlikely to be in the RIOTS4 spectral range, but the number of OB candidates is consistent with RIOTS4 having a significant incompleteness for stars in, and near, clusters.

This would also alleviate the potentially high fraction of field stars in our analysis. Our corrected observed total frequency of ejected objects is ∼36% of the entire OB population, while observations suggest field star populations closer to 20%–30%. On the other hand, the high observed ejected binary frequency (Section 2.3) supports a high ejection fraction. The reported field fractions depend on selection criteria; the youngest unbound objects, especially walkaways, are usually difficult to identify when still spatially within the cluster's projected area.

Thus, the discrepancy between the predicted and observed ratio of walkaways to runaways is mitigated by observational effects, but it does remain notable. We return to this issue below in Section 5.3.

We adopt a DES frequency of 0.33 for OB stars ejected from clusters, which yields the total observed ejection frequency of 0.36 when accounting for BSS ejections (Section 5.3). The DES runaway:walkaway branching of 30:70 adopted above yields total DES runaway and walkaway frequencies of 0.10 and 0.23 (Table 2, DES All ejected), respectively. Table 2 presents the resulting estimates for the frequencies of DES runaway and walkaway, pre- and post-SN subpopulations and binaries in the SMC, assuming the post-SN and DES binary frequencies adopted above.

Some of the 0.67 of OB stars that are not dynamically ejected from clusters are ejected by the BSS mechanism. With expected, respective frequencies of runaway and walkaway BSS ejections of 0.020 and 0.11 noted above, this yields pure BSS ejection frequencies of 0.013 and 0.074 (Table 2, Total pure BSS).

Ejected, noncompact binaries are most likely to be the tightest systems, and therefore these are destined to be progenitors of the "two-step ejection" mechanism (Pflamm-Altenburg & Kroupa 2010), which reaccelerates the surviving star upon the SN explosion of the primary. These two-step ejections are therefore a subset of both DES and BSS ejections. We note that out of our total of 11 runaway EBs and SB2s, 8 have at least one O star, increasing the likelihood that the secondary is also massive. The peak ejection age of ∼1 Myr for binaries (Oh & Kroupa 2016) is a relatively small fraction of typical OB star lifetimes (3–20 Myr), and so the vast majority of noncompact binary ejected systems will experience their first SN after the ejection event (Pflamm-Altenburg & Kroupa 2010).

The two-step process generates much higher space velocities, up to 1.5–2× faster than can be achieved by the BSS alone (Pflamm-Altenburg & Kroupa 2010). Interestingly, our fastest HMXB is the well-known object SMC X-1 (M2002-77458), which has ${v}_{\mathrm{loc}}$ = 90 ± 31 km s⁻¹ and m = 15.4 ± 1.5 M_⊙ (Figure 3; Table 1). Its PM must be confirmed, and we caution that there is inherent uncertainty in determining ${v}_{\mathrm{loc}}$ relative to other field OB stars. Paper I obtains a radial velocity for this star of 29 km s⁻¹ relative to SMC systemic, yielding a total space velocity of 95 km s⁻¹. The extreme speed relative to the small number of HMXBs is suggestive that the two-step ejection process may have played a role in this object's velocity in particular, although its measurement error remains consistent with an expected velocity around 70 km s⁻¹ for an object with its parameters (e.g., Brandt & Podsiadlowski 1995). Its short orbital period of 3.9 days (Haberl & Sturm 2016) is also consistent with a high runaway velocity.

The post-SN binaries accelerated by DES correspond to two-step ejections, and for the parameters in Table 2, their walkaway frequency is 0.017 (Table 2, DES Post-SN, binaries). Based on the DES and BSS walkaway velocity distributions (Oh & Kroupa 2016; Renzo et al. 2019), we roughly estimate that 1/3 of these two-step walkaways are reaccelerated to runaway velocities, which yields a frequency contribution 0.006 of new runaways unaccounted for by either DES or BSS models. Together with the original, runaway post-SN binaries, this gives a total two-step runaway frequency of 0.013 (Table 2, DES Two-step).

Two-step ejections may generally be observed as BSS objects. Thus, by adding the contributions from two-step runaways and pure BSS runaways, we estimate the frequency of runaways observed as BSS objects to be 0.026 (Table 2, BSS with two-step). This implies that two-step runaways may correspond to at least half of all BSS runaways. These estimates are model dependent, but in any case, we see that two-step runaways are likely a substantial fraction of BSS runaways and could dominate if the binary fraction is significantly larger than the assumed value of 30%, as appears to be the case (Sections 2.3 and 5.3). The resulting BSS population frequencies are summarized in Table 2.

With the revised values for BSS due to the effects of the two-step mechanism, Table 2 shows that the total values for predicted runaways and walkaways are only slightly adjusted to 0.10 and 0.26, respectively, from the original DES ejection frequencies (Table 2, Total predicted). The total predicted frequencies, as well as the DES/BSS ratios, are calculated using the frequencies of pre-SN DES objects and BSS objects including two-step ejections, since these correspond to what is observed. The resulting walkaway-to-runaway ratio is ∼2.6, whereas the observed ratio, corrected for walkaway incompleteness due to selection bias, is ∼1. We see that the total walkaway-to-runaway ratio is strongly dominated by the DES mechanism, and our adopted branching ratio of 70% walkaways versus 30% runaways allows about the maximum value for runaways that is plausible from the models of Oh & Kroupa (2016). Although all models predict several times more walkaways than runaways, their ratio is difficult to decrease, a problem that is exacerbated by the cluster mass function (Oh et al. 2015; Oh & Kroupa 2016). If anything, the observed ratio is too large since our 2D velocity threshold defining runaways is more stringent than the 3D threshold used in the models. As discussed in Section 5.1, accounting for possible spurious runaways and additional incompleteness due to extinction and crowding may help to resolve the discrepancy. However, our results are suggestive of a DES runaway-to-walkaway production that is higher than expected.

We also see that the total predicted noncompact binary ejection frequency is 0.21, when calculating the weighted average for both runaways and walkaways. This is low compared to the reported observed value of 0.59 for field stars (Mason et al. 2009; Lamb et al. 2016), indicating a significantly higher value for the binary ejection rate than the adopted value of 0.3. Increasing this parameter would decrease the DES-to-BSS runaway ratio, counting only pre-SN objects for DES. It would also decrease the BSS walkaway-to-runaway ratio, but these effects could be counteracted by increasing the DES runaway versus walkaway allocation, as suggested above. Increasing the total ejection frequency also slightly increases the DES-to-BSS ratio and decreases the walkaway-to-runaway ratio.

We can explore the effect of forcing the parameters to match the observations by adopting, for example, a DES runaway-to-walkaway branching ratio of 50:50 and a DES binary fraction of 1.0. This would imply that all ejected systems are expelled before their SNe or dynamical binary disruptions. We adopt a DES ejection fraction of 0.3, which yields a total ejection fraction of 0.39. Table 3 gives revised estimates for the frequencies of the various subpopulations, calculated in the same way as in Table 2. These input values produce results that are more in line with our observations: the total walkaway-to-runaway ratio has decreased to 1.2, which is much closer to the observationally derived ratio of 1.0 (Table 3), and the total ejected binary fraction is 0.58, which agrees with the observed value of 0.59 (Mason et al. 2009; Lamb et al. 2016). The DES/BSS runaway ratio is 1.8 (Table 3, DES/BSS ratio), which is still consistent with our observations, in particular, the frequency of Oe/Be stars (Sections 4.3 and 5.5).

Thus, overall, our analysis suggests a ratio of DES to BSS runaways of ∼2–3. Our data suggest that DES predictions may underestimate runaway production relative to walkaways. However, there could also be underlying issues with incompleteness and other observational biases, as described in Section 5.1. We also caution that the DES models of Oh & Kroupa (2016) and BSS models of Renzo et al. (2019) are independent, and there are likely minor inconsistencies between them, and additional physical relationships between DES and BSS ejections that are unaccounted for.

Observations of the field binary frequency more strongly suggest that the DES binary ejection fraction is higher than adopted based on the models of Oh & Kroupa (2016). This may be due to the importance of lower-mass clusters, which eject binaries at higher rates (Oh et al. 2015). Increasing the DES binary fraction also increases the frequency of two-step ejections and can do so substantially. Whereas in Table 2 the numbers of two-step and pure BSS runaways are equal, in Table 3 two-step ejections are more than 4× the number of pure BSS runaways. Thus, if the DES binary ejection fraction is indeed high, then two-step ejections likely dominate BSS runaways.

This can explain the unusually large frequency of BSS runaways (Sections 3.2 and 4.3). From fiducial models, we expect a walkaway-to-runaway ratio of ∼6 for pure BSS ejections and ∼3 including two-step ejections (Table 2), recalling that runaways are difficult to produce via BSS. But increasing both the DES binary ejection and DES runaway frequencies strongly increases the two-step contribution to the BSS runaway population, bringing the BSS walkaway-to-runaway ratio closer to the observed ratio of ∼1 (Table 3). Allowing for incompleteness in the observed number of walkaways, the contribution of two-step runaways can easily account for the high observed rate of BSS runaways.

The expected number of walkaways relative to our survey data implies that the frequency of any field OB stars that formed in situ must be small. We will address this in a forthcoming work that identifies such objects for our RIOTS4 targets by searching for associated small clusters and confirms that there are few detections (Vargas-Salazar et al. 2020). These findings support suggestions by Gvaramadze et al. (2012) that essentially all field OB stars are ejected systems.

The origin of classical Oe/Be stars has long been a puzzle and is often attributed to nonradial pulsations and/or magnetic phenomena (see, e.g., the review by Rivinius et al. 2013). However, a model that historically has received less attention is that Oe/Be stars might simply originate as objects that acquire their high rotation velocities through binary mass transfer (e.g., Pols et al. 1991; van Bever & Vanbeveren 1997; McSwain & Gies 2005; de Mink et al. 2013). Modern understanding of massive binary properties and statistics provide new leverage for this model. In particular, binary population synthesis models show that the frequency of Oe/Be runaways in the Milky Way is consistent with all of them having formed through the post-mass-transfer model (Shao & Li 2014; Boubert & Evans 2018), and observations are consistent with a prevalence of compact companions (Klement et al. 2019) and a lack of main-sequence companions (Bodensteiner et al. 2020).

We showed in Section 4.3 that the statistics and kinematics of classical Oe/Be stars in our sample are fully consistent with this population largely corresponding to post-SN secondaries that remain after the original primary has exploded. This is further supported by their statistical similarities to HMXBs. Indeed, all but one of the HMXBs in our sample are Be stars. As noted above, we find that 7 of our 14 total HMXBs, and 69 of our total 123 Oe/Be stars, have runaway velocities >30 km s⁻¹. These numbers yield walkaway-to-runaway ratios of 1.0 and 0.9, respectively, which agree well with each other, further supporting the premise that the Oe/Be stars correspond to BSS ejections. The predicted ratio of BSS walkaways to runaways estimated in Table 3 is ∼1.6, which is larger than the observed value of ∼1, a difference that is easily attributed to walkaway incompleteness, which is not accounted for. Furthermore, the predicted ratio of DES to BSS runaways in Table 3 is 1.8, which agrees with the observed ratio of 1.9 (128 DES / 69 BSS), assuming that all 69 of our runaway Oe/Be stars are from the BSS. For walkaways, the predicted ratio is ∼1.1 (Table 3), again in remarkable agreement with the observed ratio of 1.0 (53 DES/54 BSS). These statistics include the contribution of two-step ejections, which appear to dominate the BSS and, therefore, the Oe/Be population.

Renzo et al. (2019) note that most bound, post-SN systems have tight orbital periods and likely undergo mass transfer before the SN event, which thereby provides not only a simple explanation for the extreme rotational velocities but also a prediction that they must necessarily occur at relatively high frequencies. This suggests that most therefore are spun up to velocities exceeding the critical value, thus generating the excretion disks responsible for the line emission. It also suggests that most of these would also go through an HMXB phase. The consistency of the Oe/Be stars with the BSS statistics suggests that the disks are long-lived. Moreover, this model also provides a simple explanation for the strong bimodality in ${v}_{r}\sin \ i$ (Section 4.1) that is not explained by other models for the origin of Oe/Be stars (de Mink et al. 2013). Thus, while the Be phenomenon may also originate through other mechanisms, our data strongly support the post-mass-transfer model, with the vast majority corresponding to surviving, post-SN objects.