The VLT-FLAMES Tarantula Survey Observational evidence for two distinct populations of massive runaway stars in 30 Doradus

Context. The origin of massive runaway stars is an important unsolved problem in astrophysics. Two main scenarios have been proposed, namely: dynamical ejection or release from a binary at the ﬁrst core collapse. However, their relative contribution remains heavily debated. Aims. Taking advantage of two large spectroscopic campaigns towards massive stars in 30 Doradus, we aim to provide observational constraints on the properties of the O-type runaway population in the most massive active star-forming region in the Local Group. Methods. We used radial velocity measurements of the O-type star populations in 30 Doradus obtained by the VLT-FLAMES Tarantula Survey and the Tarantula Massive Binary Monitoring to identify single and binary O-type runaways. Here, we discuss the rotational properties of the detected runaways and qualitatively compare the observations with expectations of ejection scenarios. Results. We identiﬁed 23 single and one binary O-type runaway objects, most of them located outside the main star-forming regions in 30 Doradus. We ﬁnd an overabundance of rapid rotators ( v e sin i > 200kms − 1 ) among the runaway population, thus providing an explanation for the observed overabundance of rapidly rotating stars in the 30 Doradus ﬁeld. Considerations of the projected rotation rates and runaway line-of-sight velocities reveal a conspicuous absence of rapidly rotating ( v e sin i > 210kms − 1 ), fast-moving ( v los > 60kms − 1 ) runaway stars in our sample, strongly suggesting the presence of two di ﬀ erent populations of runaway stars: a population of rapidly spinning but slowly moving runaway stars and a population of fast-moving but slowly rotating ones. These are detected with a ratio close to 2:1 in our sample. Conclusions. We argue that slowly moving but rapidly spinning runaway stars result from binary ejections, while rapidly moving but slowly spinning runaways could result from dynamical ejections. Given that detection biases will more strongly impact the slow-moving runaway population, our results suggest that the binary evolution scenario dominates the current massive runaway star population in 30 Doradus.


Introduction
Massive runaway stars (Blaauw 1961) hurtle through space with speeds of tens to a couple of hundreds of km s −1 .During their lifetimes they can travel large distances (up to 1000s of pc) from their birthplaces, meaning that they can inject radiative and mechanical energy and processed chemical elements deep into the interstellar -sometimes intergalactic -medium (Ceverino & Klypin 2009;Andersson et al. 2020).
Two distinct mechanisms have been proposed to explain their origin: dynamical interaction in a dense star cluster (Blaauw 1961;Gies & Bolton 1986), or release from a close binary system when the companion explodes as a supernova (Zwicky 1957;Stone 1991).Both of these are believed to occur in nature (Hoogerwerf et al. 2000), but their relative efficiency, the resulting ejection rates and velocity distributions have been a topic of debate over the past 50 years (Blaauw 1961;Hoogerwerf et al. 2001;Jilinski et al. 2010;Dorigo Jones et al. 2020).Understanding these mechanisms and their relative contributions is pivotal to quantify the feedback of massive stars in the broader context of reionization and the chemical enrichment of galaxies (Conroy & Kratter 2012).
The 30 Doradus region (a.k.a.30 Dor) in the Large Magellanic Cloud (LMC) is the nearest strong starburst and contains one of the richest massive-star populations in the Local Group of galaxies.30 Dor therefore gives us a unique laboratory to obtain statistically meaningful observational constraints on the properties of massive stars, including those of runaways.In particular, the massive young cluster R136 (at the centre of 30 Dor) is expected to provide a suitably dense environment for dynamical interactions to occur, but is likely too young to have produce any supernova.Moreover, the metallicity of the LMC is ∼50% solar (e.g., Banerjee et al. 2012;Lebouteiller et al. 2019), so 30 Dor provides a valuable window into the properties of massive stars in conditions comparable to those during the peak of star formation in the Universe.The short lifetimes of O-type stars allow us to obtain a (close to) real-time snapshot of the production of massive runaways in 30 Dor, while also minimizing possible interlopers that would originate from outside our studied region.Here we employ measurements from multi-epoch spectroscopy of 339 massive O-type stars obtained at the Very Large Telescope (VLT) in the framework of the VLT-FLAMES Tarantula Survey (VFTS; Evans et al. 2011), combined with follow-up spectroscopy obtained over 20 months to monitor 74 of the 118 identified O-type binaries (Almeida et al. 2017).Finally, we contrast the specific line-of-sight velocities (v los ) of the identified runaways and their projected spins in an attempt to separate out their production channels.

The VLT-Flames Tarantula Survey
The bulk of the data used in this work were acquired as part of the VFTS (Evans et al. 2011), a 160h Large Programme at the European Southern Observatory (ESO).Executed from 2009 to 2010, the VFTS obtained multiepoch optical spectroscopy of 800 OB-type stars in the 30 Dor region.The targets spanned a field 20 in diameter (≈ 300 pc) centred on the R136 cluster (see Fig. 1), but avoids most of the inner 1'-region, hence most of R136 itself, due to crowding (though see Hénault-Brunet et al. 2012, for an attempt to get closer).The data were obtained with the Fibre Large Array Multi-Element Spectrograph (FLAMES, Pasquini et al. 2002), using its Medusa fibres to feed the Giraffe spectrograph with light from up to 132 targets simultaneously.The LR02 and LR03 settings of the low-resolution grating were used to provide continuous coverage of 3960-5050 Å at a spectral resolving power λ/∆λ of ∼7500.The Hα region was also observed for each target at λ/∆λ = 15 000 with the HR15N setting.
The observations and data reduction were described by Evans et al. (2011).Spectral classification, extinction properties, absolute luminosities, absolute radial velocities (RVs), projected rotational velocities (v e sin i) and atmospheric parameters of the O-type stars in the survey have been presented in various papers in the VFTS series (Sana et al. 2013;Ramírez-Agudelo et al. 2013;Maíz Apellániz et al. 2014;Sabín-Sanjulián et al. 2014;Walborn et al. 2014).Evolutionary parameters (age: τ ; current mass: M act ) were estimated for each object using a Bayesian method that matches selected observational constraints to stellar evolutionary models (Brott et al. 2011) while accounting for the observational uncertainties (Schneider et al. 2018).
O-type spectroscopic binaries and binary candidates were identified from RV variations (Sana et al. 2013).Of the 339 O-type stars observed by the VFTS, 185 showed no statistically-significant RV variations and are presumed to be single stars (though some of them will inevitably be undetected binaries; see additional considerations below on the binary detection completeness).A further 36 stars displayed significant RV variations, but with peak-topeak amplitudes of less than 20 km s −1 .These are either spectroscopic binaries or stars displaying photospheric activity (Almeida et al. 2017;Simón-Díaz et al. 2020).The remaining 118 objects, with peak-to-peak RV variations > 20 km s −1 , are considered strong spectroscopic binary candidates (Sana et al. 2013).

The Tarantula Massive Binary Monitoring
Long-term spectroscopic monitoring of 76 of the 118 Otype binary candidates was obtained over October 2012 to March 2014 in the Tarantula Massive Binary Monitoring programme (TMBM, Almeida et al. 2017).This enabled detailed characterisation of a representative fraction (65%) of the known massive spectroscopic binaries in 30 Dor.The TMBM obtained 32 spectra of each system with the LR02 setting, with the observations sampling timescales of days and weeks up to 20 months.
The data were reduced and RVs were estimated using identical methods to those employed for the VFTS data (Evans et al. 2011;Sana et al. 2013).We measured the orbital periods through a Fourier analysis and adjusted the RV curves to obtain the orbital parameters (Cumming 2004), including a refined orbital period, eccentricity, semi-amplitude of the RV curve, systemic velocity and ephemeris.Further characterisation of the double-and single-line spectroscopic binaries where reported in Mahy et al. (2020a,b) and Shenar et al. (2022a,b), respectively.
We rejected four systems with tentative orbital periods >510 days (the length of the observational campaign) as strong aliasing did not allow us to firmly derive the orbital properties.The resulting estimates of systemic velocities of the remaining 72 binaries allowed us to search for runaway binary systems using the same approach as for single stars.

Identifying runaway stars through their peculiar line-of-sight velocity
We identify runaway stars as those with estimated RVs that deviate by more than 3σ from the systemic RV of the Ostar population in the survey.Our analysis proceeds as follows.From the 339 O-type objects in the VFTS (Walborn et al. 2014), we select those without significant RV variations, i.e. the 185 presumed single stars.We first measure the systemic and 1σ dispersion of the single stars using the unweighted mean and dispersion of the RV sample and a κ − σ clipping to reject outliers.We also compute generalized RV histograms for the different O-star populations in the observed field, that we fit using Gaussian distributions (Fig. 2).The systemic velocities (v sys ) and dispersions (σ v ) obtained from these two methods and for the different O-star populations are in excellent agreement (Table 1).In doing so, we use the same definition of the NGC 2070 and NGC 2060 regions -the two main O-star associations in the 30 dor region -as adopted in the VFTS (Sana et al. 2013).These are displayed in Fig. 1.The population of O star outside the two regions are indicated as 'Remaining' in Table 1.In the following, we use the estimate of σ v obtained from the generalized histograms as they take into account individual error bars.We flag objects as runaway stars if their line-of-sight velocity (v los ) differs by more than 3σ (25.8 km s −1 ) from the systemic velocity of NGC 2070 (v 2070 sys ), which dominates the O-star population of 30 Dor (Fig. 1).We note that the mean systemic velocities of stars in the NGC 2060 and NGC 2070 differ by ∼10 km s −1 , but that all the runaways identified below are either outside NGC 2060 or have RVs  1).
that differ by more than 3σ from the systemic velocity of NGC 2060.There is thus no confusion on their runaway nature, even if NGC 2060 is their parent cluster.
The peculiar line-of-sight velocity (δv los ) of the runaway stars (both the single stars and binary systems from the TMBM) is thus given by In total, 23 (presumed) single stars and one binary meet our RV threshold (v thres = 25.8 km s −1 ), yielding an observed runaway fraction of 7 ± 1% among the O-type population of 30 Dor, where binomial statistics has been used to compute the observational uncertainties.Adopting the canonical criterion (Hoogerwerf et al. 2000(Hoogerwerf et al. , 2001) ) of a peculiar velocity larger than 30 km s −1 would reduce the number of single runaway stars by two and remove the only qualifying runaway binary system.However, it would not change the conclusions of this letter.No other runaway stars were identified in the sample with significant RV variations below 20 km s −1 .The spectroscopic binary system that qualifies as a candidate runaway system is VFTS 661 (δv los = 26.4 ± 7.1 km s −1 ).VFTS 661 is a short period (P ∼ 1.3 d) detached double-line spectroscopic and eclipsing binary.Unfortunately, the large uncertainty on its systemic velocity combined with its proximity to the RV threshold precludes a firm identification.
Most of the identified runaways have peculiar velocities of up to 60 km s −1 , but five have even larger velocities of up to 98 km s −1 (Table 2).The identified runaways are mostly located outside the two main OB associations in the region, NGC 2060 and NGC 2070 (Fig. 1).The fraction of runaway stars in these regions is 4±2% and 6±2% respectively, while it reaches 23 ± 5% in the field outside the two associations.This finding is in line with the runaway nature of the objects and strongly suggests that the massive-star field population in 30 Dor is largely composed of stars ejected from the main associations, albeit some of them with a low velocity and/or with unfavourable line-of-sight motions.The true number of runaway stars in 30 Dor is indeed expected to be much larger as three main observational effects impact our detection: the fact that our selection method is only sensitive to runaway stars that have sufficient line-ofsight velocities (v los > 25.8 km s −1 ), the finite size of our field of investigation (20 in diameter), and the completeness fraction of the survey (∼ 0.7).Contamination of our single star sample by undetected binaries would result in false detections, but this latter bias is largely outweighed by the former three effects.Detailed modelling of these biases is beyond the scope of this letter but a factor of about two to four is to be expected between the detected and true numbers of runaways in the 30 Dor region.

Projected rotational velocity distribution
We adopt the projected rotational velocities measured in Ramírez-Agudelo et al. (2013) for presumably single O stars.Typical uncertainties are of the order of 20 to 30 km s −1 .Figure 3 shows the cumulative projected spin distributions for stars identified as runaways and those that are not.A Kolmogorov-Smirnov test indicates that the two samples show differences with a significance level better than 5%.The probability of drawing by chance, from the non-runaway sample, a number of fast rotators (v e sin i > v cutoff ) as large as the one observed in the runaway sample drops well below 1% as soon as v cutoff is larger than 200 km s −1 .This again suggests that the two distributions differ significantly.The confidence of the results is consistently better than 2.5σ for v cutoff > 200 km s −1 and even reaches 3σ above 300 km s −1 .Thus, there appears to be a cutoff value above which the rotational properties of the runaway and non-runaway populations are incompatible with one another, with the runaway population being preferentially populated by rapidly spinning stars.The overabundance of fast rotators increases exponentially towards larger spinning rates.It is about a factor of two above 200 km s −1 , a factor of three above 300 and a factor of six above 400 km s −1 .Interestingly, rapid rotators are more abundant in the single O-star sample outside the NGC 2060 and 2070 OB associations than within these associations.Indeed, 32 ± 6% among the former have v e sin i > 200 km s −1 with only 20 ± 3% among the latter.It seems therefore plausible that this overabundance of fast rotators outside the main OB associations results from runaway stars spreading across the 30 Dor field-of-view.

Rapidly spinning vs. fastly moving runaways
Simultaneous considerations of the projected rotational (v e sin i) and peculiar line-of-sight (δv los ) velocities reveal distinct runaway populations (Fig. 4): a rapidly spinning slow-moving group and a slowly spinning fast-moving group.The absence of rapidly spinning fast-moving runaway stars is significant to better than 99% independent of various assumptions that can be made on the parent runaway velocity distributions.This reveals a zone of avoidance in the v e sin iδv los parameter space: a runaway desert.The existence of this runaway desert is supported by the Gaia DR3 proper motion measurements of our sample stars as discussed in the Appendix.

Discussion and conclusions
Fast rotation is one of the clearest signatures of post-binary interaction systems (Renzo et al. 2019).Under the assumption of continuous star formation, theory predicts (de Mink et al. 2013) that about 19% of the stars might be spun up to v e sin i > 200 km s −1 either through coalescence (5%) or mass-transfer (14%).A small fraction (≈ 4%) of objects are interacting but not significantly spun up (slow case-A interaction Algol-like systems).Aside from merger products (Schneider et al. 2019(Schneider et al. , 2020)), some post-interaction products may also show modest or moderate rotation, because of inclination effects or a lack of significant accretion in nonconservative mass-transfer cases.Tidal locking in relatively tight binaries may also prevent the secondary to be signifi-cantly spun up (Sen et al. 2022).The latter systems are also among those that have the largest orbital velocity so that a fraction of the high-space velocity runaways may show no significant rotational velocity increase, in qualitative agreement with the distribution of runaway stars in the v e sin i vs. δv los plane (Fig. 4).
However, the efficiency of producing fast runaways from binary evolution is model dependant (Eldridge et al. 2011;Renzo et al. 2019;Evans et al. 2020;Sen et al. 2022).Dynamical interactions occur more often in regions of high stellar density (Fujii & Portegies Zwart 2011), pointing to R136 at the core of 30 Dor as the most likely origin.The possible ongoing merger between R136 and Sabbi 1 (Sabbi et al. 2012) suggests R136 was built from subclusters with short timescales for dynamical ejection of Ostars (proportional to the relaxation timescale), possibly creating a favourable environment to further contribute to the dynamically-ejected runaway population.In either case and given its young age (Bestenlehner et al. 2020;Brands et al. 2022), it is unlikely that R136 has been ejecting stars for more than 2.5 Myr.
With 55% of the runaway O-star sample displaying v e sin i > 200 km s −1 (including the two fastest rotating stars of the entire VFTS sample), the 30 Dor runaway population presents a statistically-significant overabundance of rapidly-rotating stars compared to its non-runaway population (25% with v e sin i > 200 km s −1 ).In addition, the ejection rate of the binary channel increases with age, and the one of the dynamical channel decreases.In regards of the young age of the sample and the fact that detection biases impact comparatively more slow-moving runaways than fast-moving runaways (see Sana et al. in prep), our results support a predominance of the binary channel in producing the current and future runaway population in 30 Dor.These results contrast with the recent conclusions of Dorigo Jones et al. (2020) on the origin of runaways in the Small Magellanic Cloud.
Three of our 11 rapidly spinning -likely post binary interaction -runaway stars have masses in the range of 30 to 50 solar masses (Table 2).Ignoring dynamical effects, these high masses suggest that at least some dying stars sufficiently massive to form black holes must receive natal kicks that release their companions as runaways.Combined with the recent discovery of X-ray quiet O+BH binaries with significantly different eccentricities (Mahy et al. 2022;Shenar et al. 2022a,b), this would suggest a diversity of collapse scenarios at play.Binary-single and binary-binary exchange interactions in dense clusters can also lead to the ejection of a rapidly rotating stars, reducing the need to invoke BH kicks (Chatterjee & Tan 2012).
Finally, the overall limited contribution of dynamical interaction may result from an absence of regions with extremely high stellar density within 30 Dor.This could suggest that R136 only underwent cluster core collapse recently (or not yet), implying a moderate central density at birth.If confirmed, this would challenge the efficiency of massive star formation theories that rely on high stellar densities to significantly contribute to the formation of the most massive stars.
While the Gaia DR3 measurements support our results, the current precision does not allow us to trace back most of the runaways to their ejection locus (Platais et al. 2018), but for some rare cases (Evans et al. 2010;Lennon et al. 2018).Future Gaia data releases will hopefully provide a sufficient view of stellar motions in the three-dimensional space and provide further observational constraints on their space velocity distribution.Comparison with the observational properties of the Milky Way runaway population that is being revealed by Gaia (Maíz Apellániz et al. 2018) will provide the tools needed to investigate the impact of the environment (metallicity, parent cluster properties) on the properties of massive-star runaways.

Fig
Fig. 1.V -band Wide Field Imager (WFI) mosaic of the 30 Doradus region.O-type stars runaways single (red) and binary (green) identified in this work are spread throughout the field of view, while the non-runaway O-type stars (blue circles) cluster within the NGC 2060 and NGC 2070 OB associations (dashed 2.4'-radius circles, equiv.36 pc in radius).

Fig. 2 .
Fig. 2. Generalized line-of-sight velocity histograms of the Otype stars in 30 Dor together with the best-fit Gaussian distribution function (Table1).

Fig. 3 .
Fig. 3. Cumulative distribution of the projected rotational velocities (ve sin i) of the runaway stars (+) compared to that of the non-runaway stars (solid curve).PKS indicate the Kolmogorov-Smirnov probability that the two samples are drawn from the same parent distribution.

Fig. 4 .
Fig. 4. Diagram of the projected rotational velocities (ve sin i) versus peculiar line-of-sight velocities (δv los ) for the runaway stars identified in the 30 Doradus region.

Table 1 .
Systemic velocities and 1σ radial-velocity dispersions from the two analysis methods (see text) and runaway (RW) fractions of the single O-star populations in 30 Dor.
notes: a.At 3σ from the mean RV value of the complete and NGC 2070 samples but not for NGC 2060 (but also not spatially associated with NGC 2060).