Hunting for Runaways from the Orion Nebula Cluster

We use Gaia DR2 to hunt for runaway and walkaway stars from the Orion Nebula Cluster (ONC). We search a large area of the sky extending 45$^\circ$ in radius around the ONC and out to a distance of 1 kpc to find sources that overlapped in angular position with the cluster in the last ~10 Myr, also accounting for the ONC's own proper motion. We construct a main sample of ~17,000 runaway/walkaway candidates selected by this 2D traceback condition. However, most of these are expected to be contaminants, e.g., caused by Galactic streaming motions of stars at different distances. We thus examine six further tests or"flags"to help identify real runaways, namely: (1) possessing young stellar object (YSO) colors and magnitudes based on Gaia optical photometry; (2) having IR excess consistent with YSOs based on 2MASS and WISE photometry; (3) having a high degree of optical variability; (4) having closest approach distances to the ONC core well constrained within the half-mass radius of the cluster; (5) having direction of ejection from the ONC that avoids the main contamination zone due to Galactic streaming; and (6) having a required radial velocity (RV) for 3D overlap with the ONC that is of reasonable magnitude (or, for the small fraction, 7%, of candidate stars with Gaia measured RVs, satisfying 3D traceback). Thirteen sources, not previously noted as Orion members, pass all of these tests, while another twelve are similarly promising, except that they are located in the zone most contaminated by Galactic streaming. Among these 25 ejection candidates, there are ten with measured RVs that pass the most restrictive 3D traceback condition. We present full lists of runaway/walkaway candidate stars, estimate the high-velocity population ejected from the ONC and discuss its implications for star cluster formation theories via comparison with numerical simulations.


INTRODUCTION
The formation of star clusters takes place inside and from dense molecular clouds and is the result of gravitational collapse and fragmentation. In these early stages, young stellar clusters can be dense and the number of stars per cubic parsec increases as more stars are formed. Strong interactions between stars are expected to take place relatively frequently, especially if a significant fraction of stars are formed as binaries or higher-order multiples. Such interactions can then lead to the ejection of stars from the cluster at a range of velocities.
The Orion Nebula Cluster (ONC) is the closest massive, dense stellar cluster that is still undergoing formation, which makes it a perfect laboratory for testing star cluster formation theories. The ONC's distance is 403±7 pc (Kuhn et al. 2019). The total stellar mass is estimated to be ∼3000 M within a 3 pc radius, mixed together with an approximately similar mass of gas in this volume (Da Rio et al. 2014). However, the stars are more centrally concentrated, with their density dominating inside about 1.4 pc, having a density profile of ρ * ∝ r −2.2 extending down to 0.1 pc, where the stellar density reaches 10 4 M pc −3 (Da Rio et al. 2014).
Relatively large age spreads have been claimed to be present in the ONC (Palla et al. 2005(Palla et al. , 2007Da Rio et al. 2010, suggesting that it has been forming for the last ∼4 Myr or perhaps even longer. Three different, relatively discrete sequences in the colour-magnitude diagram have also been identified (Beccari et al. 2017), which have been interpreted as bursts of star formation interrupted by the formation and subsequent dynamical self ejection of a few massive stars (Kroupa et al. 2018). It is expected, how-arXiv:2005.04762v2 [astro-ph.GA] 12 May 2020 ever, that massive stars will also eject a number of lowermass stars. It has also been shown that the frequency and velocity distribution of ejected stars is linked to the densest state in the history of the star cluster (Oh & Kroupa 2016), as well to the dynamical timescale of its formation (Farias et al. 2019). Therefore, characterising the unbound population of a forming star cluster can provide important constraints for star cluster formation theories, including helping characterise the formation history of the system (see, e.g., Tan 2006;Tan et al. 2006).
The core of the ONC, also known as the Trapezium, is the densest and most dynamically active region of the cluster. Several runaway stars have been suggested to have been launched from this region. Using proper motions from the Hipparcos mission, the seminal work of (Hoogerwerf et al. 2001) found it highly likely that the runaway O stars µ Col and AE Aur were ejected about 2.5 Myr ago from the ONC, confirming the hypothesis first made by Blaauw & Morgan (1954). Both of these stars remain as the oldest and furthest candidate runaways to have been ejected from the ONC. High proper motion sources have also been identified in the more immediate surroundings of the Trapezium. For example, the Becklin-Neugebauer (BN) object (Becklin & Neugebauer 1967) is moving at ∼ 30 km s −1 and has been proposed to be ejected either from the θ 1 Ori C binary system in the Trapezium (Tan 2004;Chatterjee & Tan 2012), or as part of a multiple system decay involving radio sources I and a third member (Bally & Zinnecker 2005;Rodriguez et al. 2005) that is still in the Trapezium cluster region. Recent observations suggest that this other source is the nearby radio source x (Luhman 2018), which is moving at 50 km s −1 (see also Farias & Tan 2018;Bally et al. 2020).
More recently, using the high astrometric precision of Gaia, McBride & Kounkel (2019) searched an ONC membership list and identified 9 high proper motion stars. All these would have been ejected relatively recently, within the last 0.4 Myr ago. However, given that the ONC have been actively forming stars for perhaps ∼ 10 times longer, it is expected that many other high velocity runaway stars have been ejected and escaped to greater distances. However, finding such runaways becomes challenging as they mix with high proper motion field stars. Logically, most photometric efforts to identify ONC members have been focused close to the ONC. Therefore, fast runaway stars ejected more than 0.5 Myr ago and with velocities higher than 20 km s −1 that are at least 10 pc away (i.e., 1.5 • ) have likely been missed by such studies. Therefore, most of the identified ONC runaway stars are still close to the ONC region, with the notable exception of AE Aur and µ Col.
Just before submission of our paper, we have become aware of a preprint by Schoettler et al. (2020), who explored a wider area of 100 pc around the ONC, tracing back candidates using projected 2D trajectories classifying candidates by their ages which were estimated using the PARSEC isochrones from Bressan et al. (2012a). They have found 31 runaways and 54 walkaway candidates. However, they were very strict on the constraints used, i.e., by assuming an upper limit of 4 Myr for the ONC and discarding sources with 2D traceback times longer than that. They also discarded sources for which PARSEC isochronal ages were shorter than their traceback times.
In this work, we attempt to identify potential ONC runaway candidates using the unprecedented accuracy and scope of the astrometric measurements of Gaia (Lindegren et al. 2018;Gaia Collaboration et al. 2016). We explore a large area of 45 degrees around the ONC, which contains > 120 million sources in the Gaia catalog. Selecting a subset that have more accurately determined proper motions and that are within 1 kpc distance of the Sun, we trace back the two dimensional trajectories of such stars. We also test candidates under different criteria, i.e., "flags", in order to identify signatures of youth and astrometric reliability. We combine this flag system to obtain a smaller number of most probably runaway candidates. For the small fraction of the sources with radial velocity measurements, we also consider the more stringent three dimensional trace-back condition. In this way, we have produced a list of interesting sources that are prime targets for follow-up observations, e.g., to confirm if they have stellar properties consistent with ONC membership. We also discuss how the number of prime candidates that we have found compares to theoretical expectations from cluster formation simulations.

ONC FRAME OF REFERENCE
The starting point for our study is to estimate the proper motion of the ONC. A good hint of its proper motion is given by the already known runaway stars that most likely came from the ONC. Two well studied runaway stars are AE Aur and µ Col. They were first noted by Blaauw & Morgan (1954), who observed that they travel with similar velocities (∼ 100 km s −1 ) in almost opposite directions. Later, Gies & Bolton (1986) suggested that both stars were ejected from an event that involved the binary star ι Ori (see also Bagnuolo Jr. et al. 2001;Gualandris et al. 2004;Ryu et al. 2017). Hoogerwerf et al. (2001) performed a set of numerical simulations, tracing back the trajectories of AE Aur and µ Col taking into account the gravitational potential of the Galaxy. They found that the proper motion and coordinates of the star cluster that may have hosted the event are consistent with those of the ONC.
We now estimate the proper motion of the ONC from Gaia data and check how consistent it is with the results reported by Hoogerwerf et al. (2001). To do so, we make use of the membership compilation performed by Da Rio et al. . Proper motions of stars thought to be ONC members. The red circle shows our adopted proper motion for the ONC, which is obtained by taking a weighted median of the motions of ONC members. The blue star shows the proper motion of the center of mass of AE Aur, µ Col and ι Ori that are proposed to have been ejected from the ONC (Hoogerwerf et al. 2001). The green triangle shows the center of mass proper motion using the membership candidates of Da Rio et al. (2017). which also used estimated the center of mass of the ONC. We cross matched this membership list with the Gaiacatalog using the best neighbor method with a 1 cross match threshold, selecting stars within 9 from the ONC center. This angular distance corresponds to the half mass radius (1 pc) of the ONC (Da Rio et al. 2014). One of the main problems when estimating the proper motion of the ONC is the sampling. It would be tempting to use the accurate parallax measurements of Gaia to constrain the sample in the line of sight direction on an equivalent distance to the 9 used as angular distance threshold. However, if we do so the number of selected stars from the membership list is only on the order of dozens. Such proper motion would be less reliable since it would be affected by incompleteness and sparse sampling, since the number of members of the ONC is ∼ 3000 (Da Rio et al. 2014. Therefore, we have opted to use all stars flagged as members within 9 of the ONC center, which resulted in 458 member stars. The membership list compiled by Da Rio et al. (2016) also includes mass estimates of the sources, and again, it is tempting to calculate the center of mass proper motion of the selected sample. However, this quantity can be sensitive to anomalous motions of small numbers of massive stars, e.g., the most massive star system θ 1 C, which has a relatively high motion within the ONC (e.g., Tan 2004). Thus we have also measured the ONC's proper motion using a weighted median estimate. The weighted median is a robust central tendency estimator that allows the use of the individual uncertainties in proper motion. Weights are taken as 1/error 2 , and then normalized by the total sum of the weights. As in the median measurement, values are sorted by proper motion and weights are normalized by its sum. The proper motion at which the cumulative sum of the normalized weights is 0.5 is then the weighted median. Figure 1 shows the weighted median proper motion of the ONC members (red circle, hereafter adopted as the best estimate of ONC proper motion). Errors were estimated using bootstrap analysis, following the method of Kuhn et al. (2019). The resulting measurements are: µ α * = 1.43 ± 0.14 mas yr −1 and µ δ = 0.52 ± 0.12 mas yr −1 . In the figure we also show, for reference, the proper motion of the center of mass of AE Aur, µ Col and ι Ori estimated by Hoogerwerf et al. (2001) (blue star). The three proper motions are in approximate agreement.
Following the same methodology, we have measured the weighted median radial velocity of the ONC. For this purpose we have used radial velocity measurements from the INfrared Spectra of Young Nebulous Clusters (IN-SYNC) survey that covered the Orion A complex (Da Rio et al. 2017) obtaining radial velocities for 2691 sources with uncertainties in individual measurements often being 1 km s −1 . Using this catalog instead of Gaia we greatly increase the sample of sources with available radial velocities within 9 of the ONC, from 15 sources in Gaia to 200 sources marked as members of the ONC by Da Rio et al. (2016). Computing the weighted median on this sample we have obtained a radial velocity of 26.4 ± 1.6 km s −1 . This radial velocity is consistent with the one that Hoogerwerf et al. (2001) have estimated for the parent star cluster of AE Aur and µ Col of 27.6-28.3 km s −1 . In the local standard rest, this radial velocity transforms into 9.2 km s −1 . Such a radial velocity, although somewhat higher than overall average of stars in the Orion A complex of 8 km s −1 in Da Rio et al. (2017), is very consistent with that estimated from 13 CO(2-1) measurements (Nishimura et al. 2015) at the declination of the ONC (see Figure 4, 4th panel in Da Rio et al. 2017), where a radial velocity gradient is shown increasing from low to higher declinations.
Using the above estimate for the proper motion of the ONC and tracing back the trajectories of AE Aur and µ Col using great circles trajectories and assuming constant proper motions, we find that their closest approaches to the ONC are 22±28 and 44±42 respectively. These two sources are the furthest known runaway stars from the ONC. The error estimation of their closest approach comes from the errors in the proper motions. Systematic errors, such as the neglect here of the Galactic potential, will also contribute. The effects of both of these types of error grow with the traceback distance. Therefore, for sources at similar (angular) distances of µ Col and AE Aur from the ONC, we expect that true runaways will also exhibit similar errors in their closest distance to the ONC. Below, we will design our trace-back thresholds in order to capture µ Col and AE Aur.

SAMPLE SELECTION BY 2D TRACEBACK
We first select sources that are up to 45 • from the center of the ONC. Within this region there are 122,531,450 sources in the Gaia DR2 catalog. We then used the following constraints to clean the sample. First, to limit ourselves to stars with well-behaved astrometric solutions, we select those sources with reference unit weight error (RUWE) parameter < 1.4. This leaves 108,990,887 sources in the sample. Next, we select stars with parallax ( ) errors that are < 20%, i.e., /σ > 5. This reduces the sample to 18,118,187 sources. To carry out a standard variability study (described later), we require sources to have visibility periods used > 6. This makes only a minor difference to the sample size, leaving 18,113,350 sources. Finally, we restrict to sources up to 1 kpc distance (given the ONC's distance of 403 pc), which leaves a final sample of 6,760,924 sources.
After defining our sample, we wish to select new runaway candidates using variables that are common for all stars in the Gaia DR2 database. This means that we must ignore radial velocities for now, since, in our final sample, only ∼ 7% of stars have radial velocity measurements. We will do a further selection with this small subset at the end of the main analysis.
With the final sample of almost 7 million sources, we then select stars whose trajectories overlap with that of the ONC in space and time. We use the following procedure. We calculate the trajectory of each star and the ONC using assuming constant proper motion along great circle trajectories on the sky, i.e., ignoring effects of acceleration due to the Galactic potential. Each point on the trajectory has an associated trace-back time (t back ) and we use this to calculate the closest approach to the ONC (D min ) in space and time. We require that D min is smaller than a certain threshold condition given by where θ is the current angular distance of the star from the ONC. Thus the threshold becomes larger for stars that are currently further away (in angular distance) from the ONC, which allows for the fact that the errors in estimating past positions grow with longer traceback distances, i.e., due to proper motion errors, the constant proper motion approximation (which is broken by projection effects) and the effects of the Galactic potential. The normalization of the threshold condition has been adjusted to make sure that µ Col and AE Aur are recovered by this method: in particular to capture µ Col which has a closest approach D min = 43 with θ = 27 • . For sources currently close to the ONC, the threshold is ∼ 10 , which is about 1.2 pc, i.e., similar to the halfmass radius of the cluster. This choice is motivated since runaways are expected to be produced by dynamical ejection events that are more frequent in the dense, inner regions of the cluster, and having a smaller threshold helps to minimize contamination from field stars. We note that our method will select all stars that are currently within 10 of the ONC's center.
Using this trace-back method we find there are 16,994 sources that meet this 2D projected overlap condition. We considered adopting a maximum traceback time, e.g., ∼ 5 to 10 Myr, however, when we examine the distribution of traceback times of the selected sources ( Figure 2), we see that most sources are already within this range. Given the caveats of the assumed linear trajectory, constant proper motion approximation and thus possible discrepancies between 2D and 3D traceback, the real traceback time may be very different especially for sources that are far from the ONC. Thus, we simply retain all the selected sources for further analysis and note that the value of the traceback time, especially if 5 Myr, could weigh against the likelihood of a source being a genuine runaway from the ONC. Figure 3 shows the overall distribution of the ∼17,000 sources that satisfy the 2D traceback condition, described in §3. Given the number of sources detected by this method and their asymmetric distribution in position angle (PA) around the ONC, it is clear that the large majority are contaminants from Galactic field stars that have apparent past positional overlap with the ONC. It is expected that this occurs especially due to systematic differential rotation in Galactic orbits, which we will refer to as "Galactic streaming", thus explaining the asymmetry of the distribution, preferentially in  one direction from the ONC that is parallel to the Galactic plane.

SAMPLE REFINEMENT FROM SIGNATURES OF YOUTH AND KINEMATIC PROPERTIES
Given the large number of sources found by 2D traceback, the next challenge is thus to find ways to filter out most of these contaminants to identify candidates that have a higher likelihood of being real runaways. These could be targets for spectroscopic follow-up, e.g., for radial velocity measurement and better stellar characterization. To do this filtering, we now carry out six further tests, focusing on aspects of stellar youth, via: 1) optical colors; 2) IR excess colors; 3) variability; and kinematic properties, via: 4) accuracy of coincidence with the ONC centre; 5) if PA is away from the main contamination zone due to Galactic streaming; 6) radial velocity considerations to achieve 3D traceback.
If a star passes a test, we say it is "flagged" as being of greater potential interest. However, some stars are not able to be examined for all the tests, and so we will generally pass or "flag" such a star in respect of that test to err on the side of inclusion. To distinguish such cases we use a numerical value for the three states of a source with respect to a given flag: pass the flag (+1); fails (-1); and could not be tested (0).

Youth via optical color-magnitude (YSO flag)
We expect most lower-mass stars formed in the ONC to be in the pre-main sequence phase, which may be assessed from the HR diagram. In particular, we examine the position in the color-magnitude diagram (corrected by extinction and reddening) to remove contaminating lower-mass main-sequence stars, following conditions used previously by Kounkel et al. (2018); McBride & Kounkel (2019), i.e., the following cuts in colour-magnitude space in the Gaia color-magnitude system (see Fig. 4a): The extinction correction in this analysis was achieved following the method of Zari et al. (2018) for the studied region. It consists of making a 3D grid on the studied region and using the values for G band extinction, A G , and color excess E(G BP − G RP ) provided by Gaia. While the individual values for extinction and color excess are not especially ac-  curate, it is possible to use an average for sources in each bin (Zari et al. 2018;Andrae et al. 2018). We download a special sample for this purpose from the same region of our original sample, using the following conditions (Zari et al. 2018): We then grid the region in 3D using bins of 10 pc, and take the average on each bin for extinction and color excess, obtaining a 3D extinction map.
Using this method we found that 2,893 sources out of the ∼17,000 2D traceback main sample have properties consistent with YSOs (or higher-mass main sequence stars).

AllWISE IR classification (WYSO flag)
The recent study by Marton et al. (2019) performed a probabilistic classification of YSOs in the Gaia catalogue using the cross-matched table between Gaia and the AllWISE database by Marrese et al. (2019). The AllWISE source catalogue (Cutri & al. 2013) is an extension of the Widefield Infrared Survey Explorer (WISE) survey (Wright et al. 2010) that contains 747 million sources with accurate infrared photometry. WISE scanned the whole sky using four near-infrared bands at 3.4, 4.6, 12 and 23 µm, hereafter W1, W2, W3 and W4, respectively. Marton et al. (2019) used a Machine Learning approach to classify sources in four categories, i.e., Main Sequence stars (MS), Extragalactic objects (EG), Evolved stars (E) and Young Stellar Objects (YSO). For the YSO classification, they used as training sample of photometric and spectroscopic YSO catalogues listed by VizieR (see Appendix A, B and C in Marton et al. 2019, for references), Spitzer YSOs and YSO candidates from (Evans II et al. 2003). Applying a Random Forest classification, they were able to recover 93.9% of the training set correctly. Using this method they provided a probability, PL Y , of a source being a YSO using the W1-4 bands (among other features).
However, most AllWISE sources had spurious photometry and point source identification in the longest bands (W3 and W4) that could lead to false classification. Therefore the whole classification was done also using only the W1 and W2 bands, providing a complementary probability, PS Y , of being a YSO when discarding W3 and W4. Following Koenig & Leisawitz (2014) they also used random forest classification to characterize with a probability, P R , for W3 and W4 to be real. Following Unfortunately, only 4183 sources, i.e., 25%, 2D traceback sample have AllWISE photometry. Of these, we find 420 sources that fulfil the WYSO criteria (see Fig. 4b). We choose not to penalize sources that could not be evaluated using this method. This essentially means that we flag WYSO = 0 any source that was not in the Marton et al. (2019) catalog. We refer use WYSO = +1 to sources that do have AllWISE photometry and that fulfill Eq. 2.

Variability (VAR flag)
The majority of YSOs exhibit variability (e.g., Cody & Hillenbrand 2014). Thus, we expect most true ONC members will do so also. In Gaia DR2, only average values for the photometry are published for all sources, although each source has been observed at several epochs (Evans et al. 2018). The reported value of the mean flux has an associated uncertainty related to it. Variability is assessed as proportional to the standard deviation of the magnitude measurements, which can be reconstructed from the mean quantity of the flux, the flux error and the number of measurements for a given source. In Gaia, the f g band is the most precise photometric measurement, thus we use it to construct a proxy for the amplitude of the variation (Eyer et al. 2019). As amplitude proxy on g, hereafter AP g , we use its fractional standard error, i.e., AP g = σ(f g )/ f g . This value is obtained from Gaia DR2 as: where phot g n obs is the number of observations used to construct f g and phot g mean flux over error is f g divided by its error divided by √ phot g n obs. In Fig. 4c we plot AP G against f g . A noise threshold function (P A 0 ) is fit to the densest area of Figure 4c, shown with the green line, and having the form: We quantify intrinsic variability (V i ) via: This method is only well behaved for faint sources, since brighter sources suffer from other sources of photometric errors that are not well described as random noise.
Thus, we only evaluate intrinsic variability for sources with phot g mean mag = g > 13.5. Finally, we flag as VAR = +1 all sources with g > 13.5 and V i > 1. Sources with g < 13.5, which could not be evaluated using this method, are not penalized and are given a VAR = 0 flag. Overall 2,494 sources are flagged VAR ≥ 0, of which 676 are flagged VAR = 1.

Close Approach (CA flag)
Even with the unprecedented improvement in proper motion accuracy provided by Gaia, many stars still have relatively large uncertainties in their astrometric solutions. For stars that are further from the ONC, these uncertainties, plus those associated with the effects of the Galactic potential that we have not accounted for, have a correspondingly larger effect on the predicted position when the star was near the ONC. So far we have been quite generous in the closest approach distance to the ONC that is needed to select stars, i.e., it could be as large as 68.5 (∼ 8 pc) for sources that are currently 45 • away from the ONC, with this limit set to be able to recover µ Col and AE Aur (see above). Now, we wish to flag those sources that do have more accurate estimates of their proper motions that bring them within 10 of the ONC centre and with an uncertainty that no greater than 10 . Most runaways are expected to be ejected from the dense central region of the ONC, so with this flagged subset we expect to have a higher likelihood of finding real runaways and reduce the level of contamination compared to the main sample. The sources selected by this method are shown in Figure 4d: there are 1,447 sources out of the main sample of ∼ 17, 000. Thus we can see that even with this more restricted 2D traceback condition, the sample is still likely to be dominated by contaminants. We also remind that there could be true runaways, especially more distant ones like µ Col and AE Aur, that are not selected by this method.

Position Angle (PA flag)
A large degree of contamination is present after the 2D traceback selection, due to streaming in the Galactic plane. This is evident from the asymmetric distribution of the sources around the ONC (see Fig. 3). We have thus added a flag based on the position angle (PA) of a star's current angular position relative to the ONC, where 0 • is in the direction of the Galactic north pole. We flag sources that are outside the range 50 • <PA< 160 • , which is the main contamination zone. 4.6. Radial Velocity Flag (RV flag)

Stars with measured radial velocities
In the 2D traceback sample of ∼ 17, 000 sources, only about 7% have measured radial velocities, i.e., 1, 162 stars. However, for these sources, we are able to carry out a 3D . For those sources with measured radial velocities, the closest approach distance of the 3D trajectory to the ONC center, Dmin,3D, is shown versus the current 3D distance of the star from the ONC, D * ,ONC. We select those sources that satisfy the threshold condition of Eq. (8), shown by the red line. The blue symbols are AE Aur and µ Col, as labelled.
traceback analysis, which is more restrictive than the 2D method.
Using simple vector algebra, we calculate the closest approach to the ONC in 3D, D min,3D . Given the position X * and velocity V * of a star and the position X O and velocity V O of the origin, in this case the ONC, the time of the closest approach to O is Then, the closest approach distance is: Following a similar approach to that used in the 2D traceback described in §3, we adopt a closest approach threshold that grows as a function of current distance of the source from the ONC: where d ONC is the current 3D distance of a star to the ONC. The form of this equation was guided by consideration of AE Aur. While AE Aur and µ Col do not have measured radial velocities with Gaia, we made use of the values from Hipparcos used in the analysis of Hoogerwerf et al. (2001).
With these central values for radial velocity and distance for AE Aur, the past trajectory misses the ONC by 100 pc. However, we note that there are significant spreads in the distributions of the properties of the stars used by Hoogerwerf et al. (2001), while the affect of accelerations induced by the Galactic potential and Orion's giant molecular clouds' potentials will also impart apparent discrepancies.
Using the criterion described in Eq. (8) selects 516 sources as ejection candidates out of the 1200 sources with radial velocities. However, for a small number, 25, of these selected sources the obtained τ min,3D is negative, i.e., their closest approach in 3D to the ONC is in the future. We discard such sources, leaving a final selected sample of 491 out of 1,200, i.e., 40%.

Required radial velocity
In the 2D traceback selected sample 93% of sources do not have radial velocities. This is the final parameter from the 6-dimensional space needed to fully characterize the trajectory of a star. For the sources lacking radial velocity measurements, we have calculated the value of radial velocity, v r,opt needed so that their past trajectory has the closest approach to the ONC center. To account for measurement errors, we carry out Monte Carlo sampling over the distributions of astrometric parameters, assuming Gaussian distributions for the uncertainties, to obtain not only a distribution of v r,opt , but also of closest approach distances and traceback times. From these distributions we report the 16 and 84 percentiles for each source.
For sampling on , simple Gaussian sampling is not enough (Bailer-Jones 2015). Even though we have chosen to work with sources with small enough errors so we can use distances as r = 1/ , sampling distances using a Gaussian distribution around caused us to lose some sampling probability when < 0. Instead, to infer the distribution of distances that a source with a given and σ would have, we use the posterior used in Bailer- Jones et al. (2018) assuming an exponentially decreasing space density prior (see Bailer-Jones 2015;Bailer-Jones et al. 2018, for further details). Figure 6 shows distributions of v r,opt and the minimum values (16th percentile) of the distributions for D min,opt and τ 3D,opt . The first panel shows that the majority of sources need positive radial velocities to reach the ONC during their past trajectories, since most of them are at distances >410 pc. It is important to note that v r,opt generally has large uncertainties. Most of this comes from parallax uncertainties, which can give a large range of possible distances, directly affecting the radial velocity needed to reach the ONC.
The second and third panels in Figure 6 show the 16th percentile of the distributions of D min,opt and τ 3d,opt . An important point to note is that tracing back in 3D will not necessarily give the same result as 2D traceback. There are two reasons. First, in 2D traceback we assume the proper motion is constant along the trajectory, which is an approximation that becomes less valid for sources with relatively large current angular separations from the ONC. Second, the 2D traceback method does not consider the radial velocity of the ONC. For these reasons, the traceback time in 3D can be different from that in 2D.
A difference in the traceback time affects the final position in the sky of closest possible approach to the ONC. If the ONC radial velocity were zero, then the minimum closest approach from a source to the ONC would be given by the 2D closest approach in the plane of the sky. However since it moves, the closest distance may be different from its 2D counterpart. The result is that the best 3D closest approach of some sources is larger than the threshold used in the 2D traceback. We therefore exclude those sources where the minimum (16th percentile) closest approach does not satisfy the 3D traceback threshold.
Conditions on v r,opt , D min,opt and τ 3D,opt can thus be used as thresholds to exclude some sources from the candidate list. First, we consider the magnitude of v opt . We do not expect that dynamically ejected stars are likely to have radial velocities greater than 1000 km s −1 as the maximum ejection speed is approximately the escape speed from the location of ejection, which is limited by the escape speed from near the surface of the ejecting star. Indeed, known runaway stars with velocities > 100 km s −1 are very rare. If this velocity  (894) 151(2963) 269(5280) (1198) 97(278) 63 (498) 146 (637) YSO × WYSO 149 (307) 101(495) 207 (650) (111) 66 (154) YSO × WYSO × CA 50 (212) YSO × WYSO × PA 64 (138) 120 (181)  were to be used as a threshold, then 7,231 sources would be excluded. An additional 3,267 sources can be discarded if we exclude any source has min D min,opt > 2 pc, i.e., about twice the ONC half-mass radius. In principle, the revised traceback time could also be used to discard sources, however the range of values shown by τ 3D,opt are similar to those found earlier in the 2D traceback method, so we do not exclude any based on this quantity. In summary, by assessing the conditions needed for 3D traceback in the sample of sources that do not have radial velocity measurements, we have excluded 6,496 sources out of 16,994. Thus every source is given an entry for the RV flag, but with the selection criteria depending on whether it is a source with a measured radial velocity or not.

RESULTS
We now discuss the results of applying the sample refinement criteria described in the last section. Table 1 shows the number of stars that satisfy the various possible combinations of flags. For the cases of WYSO and VAR flags, not all sources could be tested, either because the sources did not appear in the AllWISE catalog (see §4.2) or were too bright (see §4.3). As mentioned in §4, we use a numerical value for each flag when passing (+1), failing (-1) or when it could not be tested (0) Table 2. Summary table of score system used to classify and label candidate runaway sources. Score is based in how many signatures of youth a source fails (indicated with -1), while the subclass column is based on how many flags for which the source could not be tested (indicated with 0).

Flag combinations
in parentheses show the combinations of flags with +1 or 0 values. In order to select the best candidates for high velocity runaways, we have developed a score system which gives an unique score to each source based on the number of flags that it meets and to the transverse velocity it has within the context of each group. There are four general scores given, as described below.
Sources that do not fail any young signature flag (i.e., YSO, WYSO and VAR) are given a score A. Sources that fail one signature of youth flag are scored with B. Sources that fail two signatures of youth are scored as C and sources that fail all three YSO, WYSO and VAR flags are scored as D. Within each score we added three subclasses depending on how many zeros are found in the signatures of youth flags. Subclass I means there is no zero in these flags, subclass II means there is one of the flags with a zero value, while subclass III means there are two flags with zero values. Therefore we have general scores as AI, AII, AIII, BI, BII, BIII, CI, CII and D.
Three modifiers are added to the score label depending on the astrometric flags CA, PA and RV. A "+" character is added if a source passes the CA flag, i.e., is a particularly good candidate whose trajectory overlaps within 10 ±10 with the ONC. A "*" character is given to sources that fail the PA flag. This reminds us of sources that are more likely to be contaminants due to Galactic streaming. Finally, a "!" character is added for sources that fail the RV flag, i.e., sources that are unlikely to come from the center of the ONC given their current astrometric parameters and uncertainties.
Finally, within each group of scores, e.g. AI+*!, an identifier is appended depending on its estimated transverse veloc-ity on the frame of reference of the ONC, v t . Where, from the fastest to the slowest candidate, a sorted ordinal number is used as identifier. Since each score is unique, we will also use it as a label for each source in the catalog. For instance, the best candidate found in our catalog is source AI+1 corresponding to the source DR2 3209590577396377856 (see Table 5). Another strong candidate, except that it happens to be in the contaminated zone, is source AI+*1 which correspond to source DR2 3015321754828860928.

Already known ONC members
Before presenting the selection of sources based on the flag system, we first discuss how the sample of known ONC members are represented in our selection system. This is motivated by the fact that from the main sample of 2D traceback selected stars, i.e., totalling 16,994sources, there are 67 A+ sources (i.e., that do not fail any flags), of which 20 are scored AI+ (i.e., pass all 6 flags), and of these samples of 67 (20) sources, 54 (17) are already known members of the Orion A complex.
McBride & Kounkel (2019) compiled a list of 5988 known Orion A complex YSO members from the literature, with 4346 of these that are part of the GaiaDR2 catalog and therefore part of the initial ∼122 million sources around 45 • of the ONC. From these, 2598 pass the clean sample criteria defined on §3 and 432 of these were found by our 2D traceback condition. Note that the 5988 sources listed by McBride & Kounkel (2019) are members of the whole Orion A complex, while the 432 sources selected by the 2D trace-back method are sources that we estimate came from the ONC or are currently within 10 of it.
We note that our 2D traceback method, in addition to finding outward moving sources, also identifies ONC members that are currently within 1.2 pc, i.e., 10 , of the ONC, but moving towards the ONC center. After correcting for perspective expansion caused by the radial motion of the ONC with respect to the Sun (see van Leeuwen 2009;Kuhn et al. 2019), the weighted median radial proper motion with respect to the ONC center for sources in this region is 0.09 ± 0.14 mas yr −1 when using the membership list of Da Rio et al. (2016), This value is smaller than that measured by Kuhn et al. (2019) (0.23 ± 0.1 mas yr −1 ). This may be due to the different membership list used, but also because we are only measuring within the half mass radius of the ONC, which should be mostly populated by bound stars. Including a wider area would also include a larger fraction of unbound stars moving outwards. Our estimated value and the distribution shown in the left panel of Figure 7 means there is no evidence for the expansion (or contraction) of the ONC center. We have checked that within this region, the ONC only shows an expansion signature when the sample is contaminated, i.e.,  (2019) with the label letters used by the authors for quick reference. Second column shows the corresponding score and label system used in this work.
Unless otherwise stated, all astrometric parameters and others derived from astrometry are taken from Gaia [J2000] (pc)  Table 3 continued rising to a maximum value of 0.24 ± 0.12 mas yr −1 if all sources in this region are used. We list the properties of highest ONC-frame velocity ONC member sources in Table 3, displaying down to 6 km s −1 , which is estimated to be twice the 2D velocity dispersion of the ONC, i.e., √ 2σ 1D , where σ 1D = 2.3 km s −1 (Da Rio et al. 2017).
The first section of Table 3 shows runaway candidates reported by McBride & Kounkel (2019). They found 9 sources coming from the center of the ONC, of which we recover 7. The two sources that we do not recover, V1916 Ori and 2MASS J05382070-0610007, were filtered out in §3 because both have high RUWE. We also recovered Brun 711 (AIII+!4), which was reported as a visitor of the ONC. Given its trajectory and very high v r,opt , we agree that it likely did not came from the ONC center. From the remaining group, 2MASS J05351295-0417499 (BI+!4), is the fastest candidate with v t = 26 km s −1 in the ONC frame. However, it does require a quite high radial velocity to reach the ONC (158 km s −1 ). The two sources that follow are AI+ sources, V1440 Ori (AI+2) and CRTS J053223.9-050523 (AI+3) which pass all flags, therefore are very likely true runaway stars. The next sources, Haro 4-379(BI+*1) and V1961 Ori (BII+9), are scored as B but are very close to passing the failed WYSO flag. Thus their chances to be true YSOs are still high, and anyway more evolved YSOs, and especially runaways, may have lost their infrared excess. Similar considerations apply to Brun 259 (BII+13), which has P(YSO) = 0.71. The remaining source V1321 Ori (BII+!29) is too bright for the variability test, although it as been classified as variable by the All Sky Automated Survey (Pojmanski 1998), but has a low P(YSO)(0.21).
From the 432 ONC members selected by the 2D traceback method, 57 have a v r,opt higher than the 100 km s −1 threshold. This only means that they likely did not come from the 10 search area used in this work, however it does not mean that they did not came from other regions of the ONC complex. The right panel of Figure 7 shows the proper motion distribution for all Orion members that were traced back to the ONC using the 2D traceback method. The source with the highest µ out is 2MASS J05430583-0807574 (BI+*!6) which moves with µ out = 152.8 ± 6.4 mas yr −1 , however it requires a very high radial velocity (-697 km s −1 ) to reach the ONC, and does not pass the RV flag criterion. In fact, many Orion members do not pass the RV flag criterion, which means that they may not have come from the ONC, but rather from other regions of the Orion complex. We also show the µ out distribution for traced back Orion members that pass the RV flag as a solid gray line of Fig. 7 right. The two fastest sources that likely came from the ONC are actually AE Aur and µ Col, that stand far from the next 3D traced

New sources
In this section we preset runaway candidates that fulfil various flags and that were not flagged as members of the ONC in the literature. Table 5 shows the best candidates from this group sectioned by score. The first three sections show sources with score AI+, AII+ and AIII+, sorted by their transverse velocity in the ONC frame within each section. The following three sections show best scored sources, but that fail the position angle flag, i.e., sources scored AI+*, AII+* and AIII+*. The trajectories of these sources are shown in Figure 8, color coded by their score.
We can see in Table 5 that the best candidate in the catalog, i.e., source AI+1, aka 2MASS J05332200-0458321, has not been flagged as being a member of the ONC in the literature. This source moves with a transverse velocity of 25.29 km s −1 in the ONC frame and was ejected ∼170,000 years ago. This source has been reported before by Rebull et al. (2000) as a M6 star, and very recently it was also reported as runaway candidate for Schoettler et al. (2020) who estimated an age of 5.0 +15 −4 Myr using PARSEC isochrones. This source is at 388 parsec form the Sun and we can see from Table 5 that its 3D distance from the ONC, d ONC , is 25 pc. While it transverse velocity is high, it is also quite close to the ONC in projection, i.e., currently at θ ONC =40'.
The other candidates that fulfil all flags and have not previously been noted as members of the ONC are DR2 3017304105587577728 (AI+7), which has a transverse velocity of 4.55 km s −1 , and DR2 3209532616814106112 (AI+17), which has v t = 1.79 km s −1 . Both sources have small transverse velocities within the recent estimates of es-cape speed of 5.6 km s −1 (Kim et al. 2019), so they are most likely still bound to the ONC.
As discussed in §4.5, the PA flag shows candidates outside the zone in the sky that is most contaminated by Galactic streaming of stars. Next, we examine sources that are positive on all criteria, but that happen to be in this contaminated zone, and therefore there is a higher chance that are false positives. There are 79 sources scored A+*, of which 32 are scored AI+*. Out of these 32 sources, only one of them is a not an already known member of the ONC. This source is DR2 3015321754828860928 (AI+*1), which moves with a high transverse velocity 80.31 km s −1 . Its current position is 3.65 degrees from the ONC and it was ejected 280,000 years ago. This source does not appear in other catalogs on Simbad, but given its luminosity is very likely to be a lowmass star. This source was also selected by the recent work of (Schoettler et al. 2020), who estimated an age of 0.4 1.1 −0.3 Myr.

Sources with measured radial velocities
About 7% of sources in our sample have measured radial velocities. None of them have been scored with the maximum score of AI+. The best source found with available radial velocity is source AII4 (2MASS J05343170-0351513), which is a known ONC YSO (Megeath et al. 2012), with a radial velocity of 20.63 km s −1 , which is within the range of the estimated v r,opt distribution to reach the ONC center, i.e., the range 19 -27 km s −1 . We note that this source is very close to passing the CA flag with D min = 7.26 ± 10.37 and therefore it likely came from the ONC.
There are three previously unknown members of the ONC scored AIII+ and 6 in the contaminated zone with score AIII+*. All these sources are too bright to be evaluated by variability and were not available on AllWISE. Source AIII+3 (TYC 4762-492-1), is the source with one of longest traceback time of the group with t back = 1.15 Myr. However, when tracing back in 3D space, we obtained a traceback time τ min,3D = 2.6 ± 0.2 Myr. Its closest approach to the ONC in 3D is 42.5 ± 9.4 pc and is currently at 186.2 pc from the ONC. Other sources with long τ min,3D are sources AIII+*7 (HD 36343, 2.4 ± 0.1 Myr) and AIII+*3 (CD-23 2974, 2.4 ± 0.03 Myr). Both of these are already classified as high proper motion stars, but with the radial velocity measured by Gaia, we estimate their trajectories approach to the ONC as close as 97 ± 6.7 pc and 103 ± 4.7 pc respectively when using straight lines trajectories. A more comprehensive analysis, including allowance for the effects of nonuniform motion, would be needed to discard or confirm such sources.
Two sources in this list are also identified by Schoettler et al. (2020) as visitors of the ONC, given their isochronal age estimates (> 20 Myr). These are sources AIII+1 and AIII+*4 (TYC 5354-1317-1). We confirm that the trajectories of these sources overlap with that of the ONC, and we estimate their closest approach distances are 0.8±6.5 pc and 15.8±10.4 pc, respectively. Source AIII+1 does not appear in other catalogs on Simbad.
6. HIGH VELOCITY DISTRIBUTION OF THE ONC 6.1. Method of estimating the observed distribution Farias et al. (2019) discussed how the statistics of the unbound stellar population created by dynamical ejection events can be an important constraint on star cluster formation models (see also Schoettler et al. 2019). However, to obtain such constraints one needs to be careful about how models are compared to observational data. While masses and ages of individual stars are challenging to determine, individual velocities rely only on astrometric measurements, such as proper motion and parallax. Unfortunately we do not have radial velocities for most of the sources from Gaia. However, with the unprecedented accuracy of parallax and proper motion measurements, we can obtain reliable transverse velocities for the ONC members and candidate runaways. In this section, we construct the high-velocity distribution of the ONC. We start this analysis, by taking the membership list used in Da Rio et al. (2016). Since such a membership list contains sources from the whole Orion complex, we select sources within 20 around the core of the ONC, which counts 740 sources flagged as members by Da Rio et al. (2016) and that are also observed by Gaia. Figure 9 shows the transverse velocity for these sources versus the different quality criteria described in §3. Most of these sources do not pass all quality criteria (shown as green areas), but the 336 sources that pass (green symbols) are part of the initial sample. From these 336 ONC members in the clean sample, 220 are captured by the 2D traceback method described in §3.
Since not all of these sources have reliable astrometry, we used the distance estimations of Bailer- Jones et al. (2018) in order to calculate transverse velocities. The resulting velocities do not change significantly for the sources that pass the quality cuts, but this step improves results considerably for sources with relatively low quality. As can be seen in Figure 9, most of the sources with v t > 10 km s −1 have large errors in their parallax (σ / > 0.2), as well as large RUWE, which means that the observations are not consistent with the astrometric model used by Gaia. Therefore, most of the sources with high v t may be caused by uncertainties in astrometry and we exclude such sources from the sample. Figure 10 shows the distribution of sources with high v t > v min . The distribution obtained when using only known ONC members that were traced-back using the method described in §3 is shown in red. There are 17 sources with velocities above 10 km s −1 in this sample with an rapid decrease in numbers as the velocity cut increases, with the fastest source V* V1175 Ori v t = 47 km s −1 .
Table 4. Sources that fullfil the RV flag with available Gaia DR2 radial velocities with scores A+. Gaia ID in bold indicate sources that are known members of the Orion A complex. v3D,ONC shows the space velocity in the frame of reference of the ONC. τmin,3D is the traceback time obtained by 3D traceback and D min,3d the corresponding 3D closest approach to the ONC.   Table 5. Gray shaded area shows the 10 search threshold around the ONC and its estimated trajectory for the past 3 Myr. With a larger shaded area showing the 20 region for reference. Symbol colors shows the different scores with AI+ (blue), AII+ (green) and AII+(red). ONC members from the literature (McBride & Kounkel 2019) are shown in black as reference. For clarity, each source is drawn once, i.e., is only shown in one of the panels. Table 5. Same as Table 3, but for sources with score A+. The table is divided in 6 sections grouped by scores AI+, AII+ and AIII+ and the same but for sources in the "contaminated" zone, i.e., AI+*, AII+* and AIII+*. A version of this table listing all ∼17,000 2D traceback candidates will be available in the electronic version of the published paper.  We have also included the corresponding distributions using the best-scored candidate runaways from this work. However, before describing such distributions, we first discuss some possible systematic uncertainties due to contaminants in our sample. As mentioned in previous sections, there is a region in the sky, which we have characterized as being from PA = 50 • to 160 • , i.e., pointing approximately parallel to the Galactic plane, where most projected Galactic orbits appear to move away from the ONC, causing an asymmetry on traceback selection and increased levels of contamination of the sample with false positives. In order to see how this contamination affects the constructed velocity distributions, we proceed with the following analysis using two approaches: Method 1: construct the velocity distribution ignoring the effects of this higher contamination zone; Method 2: construct the velocity distribution by using only the region of the sky outside this contamination zone, i.e., using only sources that pass the PA flag. Then final numbers are boosted by a statistical correction factor of 1.44 that accounts for the missing region of the sky that was not considered. During the following description, the left column of panels in Figure 10 shows results using Method 1, while the right column shows results using Method 2.
Starting with Method 1, in order to include the best candidates in the sample, we first examine all A scored candidates that pass the RV flag, i.e., with scores AI(*), AII(*), AIII(*), which includes sources that fail the PA flag. We can see from this distribution that even this sample is still likely to be highly contaminated with false positives (see top left panel in Figure 10). We consider that this distribution is unreliable for representing the true ONC high velocity population.
We next constructed another velocity distribution using only the best candidates from our sample, i.e., the sources with scores AI, but still also including sources that fail the PA flag that pass the RV criteria. This method selects a smaller fraction of candidates (black filled circles on Figure 10), which are more likely to be true runaways.
Then, the Method 1 estimate for the sample for the high velocity tail is composed by a combination of our best candidates (95 sources), 272 traced back members from Da Rio et al. (2016) that pass all the astrometric quality criteria (red filled circles on Figure 9) and we have also added AE Aur and µ Col to the final sample, that, given some overlap, has 314 sources. We show this profile as a solid line in the top left panel of Figure 10, which is the combination of the open red and black filled circles with the addition of µ Col and AE Aur that are not part of either sample.
The Method 2 estimate follows the same procedure as Method 1, described above, except excluding sources from the Galactic streaming contaminated zone. With this removal of sources that fail the PA flag, the total number of sources drops to 126. Once boosted by the correction factor of 1.44,   and compared their velocity distributions to that derived for the ONC. The selected set of models are 3000 M molecular clumps, approximated as singular polytropic spheres (McKee & Tan 2003), with a constant overall star formation efficiency of = 0.5 and different star formation efficiencies per global average free-fall time ( ff ). Models include 50% primordial binaries. Stellar evolution, including supernovae velocity kicks, is included. Simulations were performed using Nbody6++ (Aarseth 2003;Wang et al. 2015), but modified to allow the gradual assembly of star clusters including primordial binaries and a custom background potential to emulate the influence of the background gas. Depending on the model parameters, i.e., molecular cloud mass (M cl ), overall star formation efficiency, surrounding cloud mass surface density (Σ cl ) and ff , star cluster formation spans over a wide range of timescales. We have selected three sets of models that are are finished with their star formation at 1, 3 and 6 Myr. First, a low density model with Σ cl = 0.1 g cm −2 , and ff = 0.03 in which the cluster forms over 6.5 Myr. Second, a model in which the surface density is increased 10 times, Σ cl = 1.0 g cm −2 , but with the same ff = 0.03. In this case the star formation takes 1.2 Myr. For the third model, the mass surface density is also Σ cl = 1.0 g cm −2 , but ff is smaller, i.e., ff = 0.01, so that star cluster formation takes place over 3.35 Myr (see Farias et al. 2019, for further details).
Using these models we have constructed a 2D velocity distribution in order to compare with the ONC. We have constructed this 2D velocity distribution by discarding one of the velocity components in the simulations using snapshots at 1, 3 and 4 Myr. In general, we expect the high velocity population to grow as time evolves, i.e., after there has been more time for close stellar encounters leading to dynamical ejections. The bottom panel of Figure 10 shows the obtained distributions for these models. Distributions are medians of 20 realizations for each model, normalized by the current total number of member stars. The 16th and 84th percentiles are shown as the corresponding dashed lines. Note that these simulations have 50% primordial binaries and the binary fraction does not change significantly during the evolution of the cluster. For our comparison, each binary pair is counted as one single star, as it would be if such systems are not resolved, and velocities shown in the distribution are obtained from the center of mass velocities from each pair.
We compare these results with the obtained observational distributions in the previous section (solid black lines in the top panels). The total number of members of the ONC is uncertain. However for this comparison we use the total stellar mass within 3 pc estimated by Da Rio et al. (2014), of 3000 M . Using a canonical mass function (Kroupa 2001) within a range between 0.01 to 100 M , the mean stellar mass is m i ≈ 0.35 M . However, if we assume a binary fraction f bin = 0.5 and ignore higher order multiples, the mean mass per system is (see Farias et al. 2017) m s = m i (1 + f bin ) = 0.525 M . Then the number of stars including unresolved binaries in the ONC with an estimated stellar mass of 3000 M is N s = 5714. We normalize the observed distribution by the N s obtained assuming a total stellar mass between 2000 to 4500 M , i.e., with 2857 < N s < 8571. In order reach such numbers from the samples of 316 and 181 members and candidates, we complement such samples with transverse velocities drawn from a Maxwell-Boltzmann distribution with σ 1D = 2.3 km s −1 as it has been estimated in the ONC by various authors (Jones & Walker 1988;Da Rio et al. 2016;Kim et al. 2019). We show the complementary distributions used in the top panels of Figure 10, which yields the grey shaded area in the lower panels of Figure 10.
The second row panels in Figure 10 show the comparison with the low density cluster models (Σ cl = 0.1 g cm −2 ), where we can see that they are not able to produce enough high velocity runaways during such an early phases. Also, the velocity dispersion of these simulated clusters is much smaller than that of the ONC.
In the third row panels we see that these denser clusters are able to produce a similar number of high velocity stars, i.e., above 10 km s −1 , as inferred in the ONC. However, given the shorter dynamical times of these models, star formation is exhausted at 1.2 Myr. During this phase, when background gas is still present, the velocity dispersions of the simulated clusters are quite constant. However, as soon as the gas is depleted, the star cluster expands and the velocity dispersion drops, as can be seen in the low velocity end of the distribution. However, at these later times most of the strongest dynamical processing has already taken place and so the high velocity tail does not change too much after this. In the final set of models shown in the bottom row panels, star formation is less efficient and star cluster formation takes longer. While these models have a similar velocity distribution as the previous case, they are still forming stars at 3 Myr. In this respect they are a better fit to the proposed age of the ONC. However, we note that these clusters have half-mass radii of ∼ 0.3 pc at 3 to 4 Myr, which is a few times smaller than that of the ONC.
In this third set of simulations at 3 Myr the median velocity distribution is quite similar to the one estimated for the ONC, especially via Method 2, although the simulated distributions on average tend to be less populated at the highest velocities above 45 km s −1 . Still, even here the discrepancy is minor, considering the distribution in the simulations (i.e., the line of the 84th percentile) and the Poisson sampling uncertainties in the observed distribution.
Considering the global average number densities of systems inside the half-mass radius, these last set of high density simulations have a few 1000 stars pc 3 during the first ∼ 0.5 Myr, rising to ∼ 10 4 by 2.5 Myr and then declining significantly during the next few Myr, i.e., once the gas is exhausted and the cluster expands (see Figure 7 at Farias et al. 2019). Such properties may be quite similar to those expected during the formation and evolution of the ONC. However, it should be emphasized that these simulations have not been designed to match properties of the ONC. We expect within the possible parameter space of the models there will be clusters that can be a better match. The models are also relatively simple in that they so far do not assume any global elongation or spatial or kinematic subclustering of stars when they are born. Other assumptions to be investigated include effects of different initial dynamical states (currently the initial clump is in approximate virial equilibrium, including near equipartition magnetic fields), different degrees of primordial mass segregation (currently none is assumed) and the effects of primordial triples and higher order multiples (currently there are none). Overall, these observational measurements and comparisons provide a baseline to help develop and study the influence of these new ingredients in the star cluster models.

DISCUSSION & CONCLUSIONS
Using the unprecedented astrometric accuracy of Gaia DR2, we examined a wide region extending out to 45 • around the ONC to search for runaway (or walkaway) stars it may have produced. Within this area we have selected 16,994 sources that have a 2D trajectory in the sky that brings them back to the ONC in the recent past. Most of these are expected to be contaminants. Thus, using different criteria based on signatures of youth and astrometric accuracy, we have developed a scoring system that allowed us to filter and select the most likely runaway candidates. In particular, we have selected a set of 25 candidates that do not fail any of the tests, except for being in the zone of most contamination from Galactic streaming. These have not been associated with the Orion A complex or the ONC previously. Six of these sources pass all the signatures of youth tested in this work, i.e., variability, color-magnitude selection and IR YSO classification, and the fastest is escaping with a transverse velocity of 65 km s −1 in the frame of reference of the ONC. From the ∼1200 sources in the sample with available radial velocities 491 pass the most stringent radial velocity criterion to achieve 3D traceback, from which a small sample of 9 new sources (a subset of the 25 found above) do not fail any signature of youth flag, and therefore are strong runaway candidates.
Within the traceback sources, we have examined already known members of the Orion A complex. Since the trace-back method selected sources whose trajectories were consistent with the one of the ONC within its half mass radius as a threshold, current sources within this limit were all selected. We have examined the outward proper motions within this area and found no signatures radial expansion of the cluster center, with an outward median proper motion of 0.09±0.014 mas yr −1 (0.18±0.03 km s −1 ), when using only literature members in the measurement. While membership classification is somewhat challenging and a fraction of these sources may be non-member contaminants, we have shown that even if we use the whole sample, with reliable astrometric measurements, outward motions can only rise to a value of 0.23 ± 0.12 mas yr −1 (0.44 ± 0.24 km s −1 ), which is still not a clear sign of cluster expansion. This result suggests that the ONC center may have already reached dynamical equilibrium supporting the idea of a dynamically old system, even though is still in the process of forming stars (see also Da Rio et al. 2016.
Since the vast majority of sources in our traceback sample do not have measured radial velocities, we have computed the distributions of radial velocities that minimize the closest 3D trajectory to the ONC. We have obtained these distributions by sampling each astrometric quantity within their errors using a Monte Carlo approach. This derived sub-product of this work can be very useful for quickly distinguishing good candidate runaways when new RV measurements are obtained in the future. However, even though we have only used sources with the best astrometry, there are still a considerable group of sources whose range of optimal velocities are extremely large, and for which this metric does not provide very strong constraints on an ejection scenario. The selection criterion we used with this quantity was then rather restrictive, since we marked as negative sources with required radial velocity above 100 km s −1 . Such restrictive threshold was designed to eliminate most false positives, however some runaway stars with larger velocities are still possible and may still be hidden in our sample.
We have also estimated the total high velocity distribution of the ONC using the known and new members with the best astrometry and membership probability, selecting a sample of about 200 to 300 sources, depending on the method. While there are still significant systematic uncertainties in the estimation of this distribution, we have compared it with theoretical models based on simulations that include realistic fractions of primordial binaries and gradual formation of stars, which are necessary components for accurately capturing rates of dynamical ejections. These simulations can successfully reproduce the normalization and shape of the estimated velocity distribution of the ONC, however only when using higher density models (Sigma c l = 1g cm −2 ) with relatively slow star formation ( ff = 0.01). A more general exploration of simulation parameter space for clusters specif-ically designed to match the ONC is the next step to be conducted.
Very recently we have become aware of a parallel work from (Schoettler et al. 2020), who have performed a similar observational search for runaway stars from the ONC, but restricting to a smaller region extending only 100 pc away, i.e., up to 14 • from the cluster. They also performed N-body simulations of clusters, which in their case were initialised with sub-virial, fractal distributions of stars with a primordial binary fraction that depends on primary mass (∼50% on average). While the ultimate goal of their study is similar to this work, their methods of selection and classification of sources are very different. Some of the main differences are that they base the classification in the comparison with the traceback time and the estimated ages for their sources using isochrone fitting from PARSEC (Bressan et al. 2012b). They were quite restrictive on the age of the ONC, excluding from the runaway lists sources with ages larger than 4 Myr. While this is a reasonable limit for most ONC stars (Da Rio et al. 2016), we consider that it may be too restrictive a choice given the uncertainties in age estimates of individual stars from isochrone fitting. Schoettler et al. (2020) have reported 54 walkaway candidates with velocities between 10-30 km s −1 and 31 runaway candidates with velocities above 30 km s −1 . From their sources, we recovered 22 out of 31 runaways and 32 out of their 54 walkaways, but give them varying degrees of likelihood of being true runaways. On the other side, from our 25 best candidates they have listed 4 as being ejected from the ONC and 2 as visitors. These differences in the number of selected sources, given that we observed a wider area, likely arise from the fact that we do a more restrictive initial astrometric cleaning and that we also use a more restrictive baseline boundary for the traceback method, i.e., the half mass radius of the ONC (1.2 pc) compared to the 2.5 pc boundary used in their work. Two of their walkaway and two of their runaway stars are among our top ranked candidates. These walkaways are sources AI+1, which is our fastest source that passes all the flags (their estimated age for this source is 5 15 −4 Myr), and AIII+4 that we could not evaluate on variability or YSO probability (their estimated age for this source is 0.8 9 −0.7 Myr). The coinciding runaway candidates are sources AI+*1 and AII+*1, that fulfill all quality criteria, but that we note are in a zone contaminated by Galactic streaming. Based on age, they have flagged as candidates two sources that we confirm likely come from the ONC based on 3D traceback, these are sources AIII+1 and AIII+*4, for which they have estimated ages larger than 20 Myr. The majority of the sources presented here, together with most of the candidates presented by Schoettler et al. (2020), have missing radial velocities. One of the most important follow-up observations stimulated by this work thus involves radial velocity measurements of our candidate runaways to confirm their ejection with the more restrictive 3D traceback criterion.
We remark that one important goal of finding the oldest ejected runaways, which may be currently hidden among a cloud of Galactic field contaminants, is to constrain the star formation history of clusters, i.e., the timescale of star cluster formation. For the ONC, this runaway age constraint is still set by µ Coland AE Aur, ejected about 2.5 Myr ago (Hoogerwerf et al. 2001). The formation time is obviously a basic parameter for cluster formation theories (e.g-. Tan 2006;Nakamura & Li 2007) and also influences estimates of fundamental star formation properties, such as the efficiency per freefall time (e.g., Da Rio et al. 2014). The cluster formation time is also directly related to the formation timescale of massive star formation in competitive accretion models (Wang et al. 2010) as discussed by Tan et al. (2014). Thus a dedicated search to find relatively old, likely lower-mass runaways from the ONC should be attempted, with a starting point from our main 2D traceback sample.