Stellar Proper Motions in the Orion Nebula Cluster

The observations from HST consisted of 11 epochs between 1995 and 2015 (Prosser et al. 1994; O'dell et al. 1997; Rubin et al. 1997; O'Dell 2001; Robberto et al. 2004, 2013), mostly with medium or wide optical/IR passband filters (F435W, F439W, F539M, F555W, F775W, F791W, and F139M), except for IR filter F130N. All HST archival images having central coordinates within ∼30 of the center of the ONC were selected. However, only those images with exposure times longer than 40 s were used in our PM analysis to ensure sufficiently high signal-to-noise ratios for faint sources. A few of these images were also rejected in the process of matching and alignment (see Section 3.2).

The HST images were obtained from several different cameras. The ACS/WFC consists of two 2048 × 4096 pixel CCD detectors. The plate scale is 50 mas pixel^–1, which corresponds to a 202'' × 200'' field of view. The WFPC2 uses four 800 × 800 pixel CCDs where three of them cover a 150'' × 150'' region (WF) and have a pixel scale of 100 mas pixel⁻¹. The fourth CCD (PC) images a 34'' × 34'' field with a spatial scale of 56 mas pixel⁻¹. The WFC3IR channel uses a single 1024 × 1024 pixel CCD detector with a plate scale of 130 mas pixel⁻¹, corresponding to a 136'' × 123'' field of view.

Observations that were within ∼1 month and with the same instrument were combined to define a single epoch. In Table 1, we provide the complete list of HST observations for the different epochs used in this work, including the epoch number, dates of observations, R.A. and decl. at the center of the frames, instrument, filter, total exposure time, and principal investigator for the data.

The observations with NIRC2 (instrument PI: K. Matthews) focused on the BN/KL star-forming region (α = 05:35:14.16, δ = −05:22:21.5). The data were obtained on 2010 October 30–November 1 and 2014 December 11–12. The first run in 2010 is also described in Sitarski et al. (2013). The observations were conducted using the laser guide star adaptive optics (LGS-AO) system (Wizinowich et al. 2006). The LGS corrected for most atmospheric aberrations; however, low-order tip-tilt terms were corrected using visible-light observations of the star Paranego 1839 (${\alpha }_{{\rm{J}}2000}$ = 05:35:14:64, δ_J2000 = −05:22:33.7). In order to avoid the strong nebulosity in this region, sky frames were obtained for the wavefront sensors using larger-than-normal sky offset positions.

The two epochs of Keck AO observations covered nearly the same sky area with the wide-field cameras on NIRC2, which has a pixel scale of 39.686 mas pixel⁻¹ (Yelda et al. 2010) and field of view of 40'' × 40'', in the same passband, He i b (λ₀ = 2.06 μm, Δλ = 0.03 μm). The narrowband filter allows us to avoid the saturation of bright sources such as BN. The images were mosaicked around the BN/KL region for a total areal coverage of 1.4 arcmin². Sky frames were taken interspersed with science observations in a dark region ∼15° to the east. Sky observations were timed in such a way that the field rotator mirror angle was identical to that of the science exposures, which is necessary to accurately subtract thermal emission from the field rotator mirror in this band (Stolte et al. 2008).

A summary of our Keck AO observations is listed in Table 2. The field of view of our Keck data is illustrated by a dashed polygon in Figure 1.

Zoom In Zoom Out Reset image size

3.1.1. HST

3.1.2. Keck AO

The first step in measuring relative PMs is to establish an astrometric reference frame. With distortion-free astrometric measurements, Gaia DR2 (Gaia Collaboration et al. 2016, 2018) provides an excellent foundation for a reference frame. One caveat of using Gaia DR2 alone as a reference frame is that in the central region of the ONC, the photometry only reaches G ∼ 17, mostly due to high nebulosity. Given the star catalogs from different passbands, we needed a reference frame for alignment that covers a wider range of magnitudes and colors.

We constructed a new reference frame, building upon Gaia DR2 using our star catalogs from ACS/WFC F775W and WFC3IR F139M in epochs 8–11 as follows. First, we converted the celestial coordinates of stars in Gaia DR2 into right-handed Cartesian coordinates x and y parallel to the R.A. and decl. directions (i.e., $x={\rm{\Delta }}\alpha \cos \delta$, y = Δδ),¹⁰ respectively. In order to compromise between the number of reference sources and the accuracy of the astrometry, we adopted only Gaia stars with measurement errors smaller than 60 mas. We cross-matched the stars in the HST F775W and F139M catalogs with those of the Gaia DR2 catalog within a radius of 60 mas and applied 2D second-order polynomials to transform the positions of stars in the input star catalogs to those of the reference stars. Although the minimum number of data points for a second-order polynomial fit is six, we excluded frames that had less than nine matches to ensure a good fit. The median position was calculated for each star from both Gaia DR2 and the transformed catalogs. During each polynomial fitting, we rejected 3σ outliers and repeated the same process until convergence into a final solution.

In the final iteration, we used the new reference frame to transform and match the rest of the HST and Keck star catalogs from all epochs. A final catalog was constructed using the median of the transformed positions for each star in each epoch, and the standard error of the median was adopted as the positional uncertainty.

The relative PMs of each star were measured using least-squares straight-line fits to the positions over time. We set a lower limit for the time baseline (Δt) as 1 yr. For stars identified in three epochs or more, the linear fit determined the velocities with the errors calculated using the covariance matrix from the fit. The velocities of the remaining objects, detected in only two epochs, were calculated as the positional difference between the two epochs divided by the time baseline, and their uncertainties were calculated as the quadratic sum of positional errors in each epoch. We provide our PM catalog with positions at epoch 2015.5 in Table 3. To estimate the photometric depth of the catalog, we cross-matched the PM catalog with the first-pass hst1pass outputs for F139M images, lowering the flux limit to 20 electrons. We found that our sample includes 693 stars that fall within the range between F139M = 9.5 and 20.5, ∼95% of which are brighter than F139M = 18.0, and 22 stars that are undetected in the near-IR passband.

Table 3.  PM Measurements

ID	αJ2000a	δJ2000a	${\mu }_{{\alpha }^{* }}$	${\epsilon }_{{\mu }_{{\alpha }^{* }}}$	μδ	${\epsilon }_{{\mu }_{\delta }}$	Ndet	Δt	F139	Noteb
	deg	deg	mas yr−1	mas yr−1	mas yr−1	mas yr−1		yr	mag
1	83.85346518	−5.36602620	0.34	0.09	−0.71	0.05	4	16.9	14.57
2	83.83552846	−5.34779414	−0.20	0.16	0.68	0.55	3	16.3	12.12
3	83.77062885	−5.35257430	−0.43	0.06	2.12	0.18	3	11.0	13.25
4	83.77579375	−5.33816393	−0.41	0.28	−0.75	0.30	3	11.0	14.52
5	83.82285598	−5.38091090	1.38	0.16	1.64	0.05	5	16.4	13.35
6	83.79410989	−5.37909779	0.14	0.46	−1.53	0.06	3	16.3	12.11
7	83.81797732	−5.34033710	−1.00	0.13	0.01	0.45	3	16.3	12.93
8	83.85596135	−5.39258964	−0.79	0.11	−0.44	0.11	4	16.3	12.81
9	83.84994021	−5.41941865	0.12	0.81	0.98	0.15	4	14.6	13.04
10	83.79543932	−5.37954202	−0.04	0.27	1.00	0.40	5	19.4	13.28
⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯	⋯

Notes.

^aEpoch 2015.5. ^bIn this column, 1 = proplyd morphology; 2 = Herbig–Haro object (Ricci et al. 2008); 3 = double star, previously known in the literature (Hillenbrand 1997; Robberto et al. 2013); 4 = new double star, previously reported as single in the literature; 5 = double-star candidate, previously reported as single in the literature; 6 = double-star candidate, previously unreported.

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

Download table as:  DataTypeset image

From the measured PM measurements, we derived the internal PM dispersion of the ONC using Bayes's theorem and a multivariate normal distribution model. Assuming that each PM measurement is drawn from a Gaussian distribution with mean $\bar{\mu }$ and an intrinsic dispersion σ_μ, the likelihood for the ith PM measurement μ_i ± _i is then given by

$$\begin{eqnarray}&&{L}_{i}({\mu }_{i}| \bar{\mu },{\sigma }_{\mu })=G({\mu }_{i}| \bar{\mu },\sqrt{{\sigma }_{\mu }^{2}+{\epsilon }_{i}^{2}}),\end{eqnarray} \tag{ 1 }$$

where the final dispersion is the quadratic sum of the intrinsic dispersion, ${\sigma }_{\mu }$, and the error on the measurement for each star, _i. Given a set of measurements $D\equiv \{{\mu }_{i},{\epsilon }_{i}{\}}_{i=1}^{N}$, the posterior probability P is defined by Bayes's theorem as

$$\begin{eqnarray}&&P(\bar{\mu },{\sigma }_{\mu }| D)\propto \displaystyle \prod _{i}{L}_{i}({\mu }_{i}| \bar{\mu },{\sigma }_{\mu })p(\bar{\mu },{\sigma }_{\mu }),\end{eqnarray} \tag{ 2 }$$

where $p(\bar{\mu },{\sigma }_{\mu })\equiv p(\bar{\mu })p({\sigma }_{\mu })$ is the prior for the mean and the standard deviation. We adopt a flat prior for the mean and a "noninformative" Jeffreys prior for the standard deviation, i.e., $p({\sigma }_{\mu })\propto {\sigma }_{\mu }^{-1}$ (see Section 7 in Jaynes 1968 for justification).¹¹ As each star has PM measurements along R.A. (α) and decl.(δ), we maximized the product of the posteriori for the two axes, i.e., ${P}_{\alpha ,\delta }={P}_{\alpha }{P}_{\delta }$, utilizing the Markov chain Monte Carlo (MCMC) Ensemble sampler emcee (Foreman-Mackey et al. 2013).

We compare our PM measurements with those in Gaia DR2. Due to the strong nebulosity, the astrometric measurements from Gaia DR2 have overall lower quality for stars around the ONC compared with other nebula-free regions. As a result, we adopt a generous quality cut for Gaia DR2 stars in order to compromise between the astrometric quality and the size of the comparison sample¹² : astrometric_gof_al<16 and photometric_mean_g_mag<16. With this condition, we found 15 matches between Gaia DR2 and our PM catalog. The main panel of Figure 3 shows the difference in PMs along the R.A. and decl. axes, where the data points are concentrated around (0, 0) within 1 mas yr⁻¹. Figure 4 verifies that the PM vectors tend to point in similar directions with similar magnitudes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Returning to the main panel of Figure 3, however, we noticed that the differences generally exceed the measurement errors and exhibit an asymmetric distribution. The inconsistency in the amplitudes of PMs can be attributed to underestimated PM uncertainties in Gaia DR2. Arenou et al. (2018) demonstrated that the parallax and PM errors in Gaia DR2 are underestimated and tend to overestimate the intrinsic dispersions for distant open/globular clusters. To test the possibility of underestimated uncertainties in Gaia DR2 around the ONC, we compared the PM dispersions derived from the stars in common between Gaia DR2 and our catalog, excluding kinematic outliers identified in Section 4.2. First, we compared the PM dispersions of the 15 stars in Figure 3. The Gaia DR2 resulted in PM dispersions (${\sigma }_{{\alpha }^{* }},{\sigma }_{\delta }$) = (${1.33}_{-0.23}^{+0.31},{0.91}_{-0.15}^{+0.21}$) mas yr⁻¹, which is nearly 30% larger than the dispersion from our PM measurements, (${\sigma }_{{\alpha }^{* }},{\sigma }_{\delta }$) = (${1.04}_{-0.19}^{+0.26},{0.72}_{-0.13}^{+0.17}$) mas yr⁻¹.

In order to reach comparable PM dispersions, the PM uncertainties in Gaia DR2 need to be increased by a factor of ∼3. Given the small size of the sample, we performed the same test using a larger sample, without the quality cut we applied for the previous sample. We found 197 matches between Gaia DR2 five-parameter sources and our catalog after excluding the kinematic outliers listed in Section 4.2 and measured the PM dispersions following the same procedure described in Section 4.3. We find a 2D dispersion from the Gaia DR2 measurements of (${\sigma }_{{\alpha }^{* }},{\sigma }_{\delta }$) = (1.50 ± 0.09, 1.69 ± 0.10) mas yr⁻¹, which is nearly 70% larger than the dispersion for our catalog of (0.83 ± 0.05, 1.00 ± 0.06) mas yr⁻¹. For comparable PM dispersions, increasing the PM uncertainties in Gaia by a factor of ∼3 is again required. We note that this increased error is still needed even when comparing to the PM dispersions for our entire sample or from previous surveys at optical wavelengths (see Section 4.3 and Jones & Walker 1988). The inset of Figure 3 is the same as the main panel but with PM uncertainties in Gaia DR2 increased by a factor of 3, where the differences appear to be consistent overall with zero within ∼1σ. There is one outlier whose PM difference is offset toward the southeast. This star has the smallest number of both good observations and visibility periods in Gaia DR2 out of all the matched stars, whose astrometric_n_good_obs_al and visibility_preiods_used values range from 106 to 205 and from 7 to 12, respectively. The outlier thus falls within the regime of possible systematic errors in PMs of Gaia DR2 induced by the scanning law of the survey as demonstrated in Appendix A. The direction of the bias is also consistent with the scan direction around the ONC, northwest–southeast, which can be traced by the positions of Gaia DR2 sources filtered based on the number of observations or visibility periods. We therefore conclude that, compared to Gaia DR2, there is no significant discrepancy ascribed to our PM measurements.

Our analysis recovered known high-velocity stars, BN and source x, in the BN/KL complex. We note that these sources were not detected in the optical but only in the IR images from Keck/NIRC2 and WFC3/IR. We measure PMs of (${\mu }_{{\alpha }^{* }}$, μ_δ) = (−7.2 ± 2.7, 12.2 ± 1.9 mas yr⁻¹) for BN and (26.8 ± 1.5, −18.4 ± 1.5 mas yr⁻¹) for source x. Our measurements agree with those both at IR wavelengths from Luhman et al. (2017) and in the radio from Rodríguez et al. (2017) within 1σ.

We also recovered source n in the BN/KL complex, whose PM was reported to be as high as ∼7 mas yr⁻¹ in some radio studies (Rodríguez et al. 2017). In fact, the source appears highly elongated at radio wavelengths, which hinders a reliable PM measurement. At IR wavelengths, it appears as a single point source with a much smaller PM value (Luhman et al. 2017). Not surprisingly, our PM measurement for source n (1.9 ± 1.0, 1.0 ± 0.7) mas yr⁻¹ agrees reasonably well with the motions of (−1.8 ± 1.4, −2.5 ± 1.4) mas yr⁻¹ previously measured in the IR by Luhman et al. (2017) or (1.6 ± 1.6, 3.4 ± 1.6) mas yr⁻¹ in the millimeter (Goddi et al. 2011) but disagrees with the PM of (0.0 ± 0.9, −7.8 ± 0.6) mas yr⁻¹ in the radio data (Rodríguez et al. 2017).

In addition to the previously known high-velocity stars, we detected three other stars with large PMs, as shown in the left panel of Figure 2. We note that we initially had a few more candidate stars with large PM values; but, after visual inspection, they were identified as false positives ascribed to marginally resolved double stars, Herbig–Haro objects, or proplyds (Prosser et al. 1994; Hillenbrand 1997; Reipurth et al. 2007; Ricci et al. 2008; Robberto et al. 2013; Duchêne et al. 2018). Since foreground field stars often have large PMs, we performed further investigation to determine the nature of the three objects. Figure 5 shows the HST/ACS photometry of stars covering ∼600 arcmin² around the ONC from Robberto et al. (2013). The colors and magnitudes of the four objects are systematically bluer than the young ONC sequence and found in the locus of foreground objects with low reddening. This finding strengthens the idea that these kinematic outliers are most likely field stars. We note that the brightest object, shown by a cyan pentagon in Figure 5, corresponds to source 583 in the photometric survey of Hillenbrand (1997), where the star was classified as a nonmember with 0% membership probability. This object also corresponds to source 3017360902234836608 in Gaia DR2 with a relatively large parallax, 4.14 ± 0.07 mas (≡242 ± 4 pc), which supports that it is likely a foreground star. For internal kinematic analysis, we do not include these three peculiar objects. Also excluded are BN and source x, as was done in previous studies (e.g., Jones & Walker 1988; Dzib et al. 2017).

Zoom In Zoom Out Reset image size

We also identify probable escaping, or evaporating, stars whose high velocities deviate significantly from our Gaussian velocity distribution model (see Section 5.1). As demonstrated by Kuhn et al. (2019), the deviation can be visualized on a plot of observed data quantiles (Q_data) versus the theoretical quantiles of the Gaussian distribution (Q_theo). The quantiles are

$$\begin{eqnarray}&&{Q}_{\mathrm{data},i}=\displaystyle \frac{{\mu }_{i}-\bar{\mu }}{\sqrt{{\sigma }_{\mu }^{2}+{\epsilon }_{i}^{2}}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{Q}_{\mathrm{theo},i}=\sqrt{2}\,{\mathrm{erf}}^{-1}(2({r}_{i}-0.5)/n-1),\end{eqnarray} \tag{ 4 }$$

where ${\mathrm{erf}}^{-1}$ is the inverse of the error function, n is the number of measurements, and r_i is the rank of the ith measurement. The mean $\bar{\mu }$ and standard deviation σ_μ are computed with the method described in Section 3.3. The upper panels in Figure 6 show the Q–Q plots of all stars except for the high-velocity stars. Overall, the velocity distribution is well fit by a normal distribution along both the α and δ axes. In the α axis, however, we notice that beyond ${Q}_{\mathrm{data},{\alpha }^{* }}=\pm 3$, the data quantiles deviate from the expected quantile function of the Gaussian distribution. In the lower panels, after excluding the six stars outside ${Q}_{\mathrm{data},{\alpha }^{* }}=\pm 3$ (red filled circles) and recalculating the mean and the standard deviation, the velocity distribution exhibits consistency with a normal distribution within the 95% confidence envelope.

Zoom In Zoom Out Reset image size

We consider the possibility that the outlier stars are evaporating from the cluster by comparing their velocities to the escape speed. Using the virial theorem, the mean-square escape speed can be estimated as

$$\begin{eqnarray}&&\langle {v}_{e}^{2}{\rangle }^{1/2}=2\langle {v}^{2}{\rangle }^{1/2},\end{eqnarray} \tag{ 5 }$$

where $\langle {v}^{2}{\rangle }^{1/2}$ is the mean-square speed of the cluster's stars (Binney & Tremaine 2008). Using this relation, we approximate a corresponding angular escape speed of ≈3.1 mas yr⁻¹. We found that the apparent angular speed of the outliers ranges from 3.2 to 4.5 mas yr⁻¹, all of which exceed the speed limit for evaporation. We identified 12 additional candidate stars from our sample whose apparent angular speeds exceed this threshold (see Figure 2), although they do not stand out as statistically significant outliers on the Q–Q plots partially due to large measurement errors or dispersion along the δ axis (see Figure 2). We note that we initially had several more candidates that were excluded after visual inspection showed that they were marginally resolved double stars or unresolved double-star candidates with highly elongated, double-lobed morphology in HST/ACS or Keck/NIRC2 images. Hereafter, we refer to the relatively high- and low-significance escaping star candidates (ESCs) as ESC group 1 and ESC group 2, respectively. In Section 5.1, we demonstrate that these stars preferentially occupy the central region of the ONC. To account for their effect on the radial variation of the velocity dispersion, we present the PM dispersion as a function of radius in Table 5 for three cases: excluding (a) none of ESC groups, (b) ESC group 1, and (c) ESC group 1 + 2. Otherwise, we only exclude ESC group 1 when modeling the PM distribution in the following sections.

We note that all false positives are included and flagged in our PM catalog (Table 3). Among the false positives, we have newly identified two double stars and two double-star candidates.

Using the method described in Section 3.3, we obtain the mean PM and intrinsic dispersion of the ONC in Cartesian coordinates:

$$\begin{eqnarray*}\begin{array}{rcl}{\bar{\mu }}_{{\alpha }^{* }} & = & -0.04\pm 0.03\,\mathrm{mas}\,{\mathrm{yr}}^{-1},\\ {\bar{\mu }}_{\delta } & = & -0.05\pm 0.05\,\mathrm{mas}\,{\mathrm{yr}}^{-1},\\ {\sigma }_{\mu ,{\alpha }^{* }} & = & 0.83\pm 0.02\,\mathrm{mas}\,{\mathrm{yr}}^{-1},\\ {\sigma }_{\mu ,\delta } & = & 1.12\pm 0.03\,\mathrm{mas}\,{\mathrm{yr}}^{-1}.\end{array}\end{eqnarray*}$$

Following the same procedure, we also computed the mean PM and dispersion along the radial axis away from the cluster center and the tangential axis perpendicular to it:

$$\begin{eqnarray*}\begin{array}{rcl}{\bar{\mu }}_{r} & = & -0.04\pm 0.04\,\mathrm{mas}\,{\mathrm{yr}}^{-1},\\ {\bar{\mu }}_{t} & = & 0.00\pm 0.04\,\mathrm{mas}\,{\mathrm{yr}}^{-1},\\ {\sigma }_{\mu ,r} & = & 0.97\pm 0.03\,\mathrm{mas}\,{\mathrm{yr}}^{-1},\\ {\sigma }_{\mu ,t} & = & 1.00\pm 0.03\,\mathrm{mas}\,{\mathrm{yr}}^{-1}.\end{array}\end{eqnarray*}$$

Our PM dispersions agree with those found by Jones & Walker (1988): (${\sigma }_{\mu ,{\alpha }^{* }},{\sigma }_{\mu ,\delta }$) = (0.91 ± 0.06, 1.18 ± 0.05), (${\sigma }_{\mu ,r},{\sigma }_{\mu ,t}$) = (1.06 ± 0.05, 1.04 ± 0.05) mas yr⁻¹.

The PM dispersions were also measured for stars grouped into equally partitioned magnitude bins in F139M (N = 110) and by distance from the center of the ONC, given in Tables 4 and 5. The one-dimensional PM dispersions ${\sigma }_{\mu ,1{\rm{D}}}$ were obtained by taking the quadratic mean of R.A. and decl. dispersions, ${\sigma }_{\mu ,{\alpha }^{* }}$ and ${\sigma }_{\mu ,\delta }$. The PM dispersions appear to be essentially flat within the uncertainties from ${m}_{{\rm{F}}139{\rm{M}}}=9.50$ to 16.09, below which there is marginal evidence of decreasing R.A. dispersion. The PM dispersions more obviously decrease with radius from the center to $R=3\buildrel{\,\prime}\over{.} 0$.

Table 4.  PM Dispersions as a Function of Magnitude

F139M	${\sigma }_{\mu ,{\alpha }^{* }}$	σμ,δ	${\sigma }_{\mu ,r}$	σμ,t	${\sigma }_{\mu ,1{\rm{D}}}$	N
mag	mas yr−1	mas yr−1	mas yr−1	mas yr−1	mas yr−1
9.50–12.18	0.86 ± 0.07	1.08 ± 0.08	0.91 ± 0.07	1.02 ± 0.08	0.98 ± 0.05	110
12.18–13.00	0.87 ± 0.07	1.18 ± 0.09	1.09 ± 0.08	1.00 ± 0.07	1.04 ± 0.06	110
13.00–13.56	0.84 ± 0.07	1.11 ± 0.09	0.92 ± 0.07	1.05 ± 0.08	0.98 ± 0.06	110
13.56–14.45	0.85 ± 0.07	1.19 ± 0.09	1.01 ± 0.08	1.07 ± 0.08	1.03 ± 0.06	110
14.45–16.09	0.83 ± 0.06	1.12 ± 0.08	1.06 ± 0.08	0.94 ± 0.07	0.99 ± 0.05	110
16.09–18.38	0.67 ± 0.05	1.03 ± 0.08	0.79 ± 0.06	0.94 ± 0.07	0.87 ± 0.05	110

Note.

^aThe dispersion columns give 1D intrinsic dispersions in the R.A., decl., radial, and tangential directions. The final dispersion column is the mean of the R.A. and decl. dispersions.

Download table as:  ASCIITypeset image

Kuhn et al. (2019) demonstrated that the ONC is one of only a few young open clusters whose stellar velocities are consistent with a multivariate normal distribution or a thermodynamic Maxwell–Boltzmann distribution (see Figure 10 in Kuhn et al. 2019). For this analysis, however, the authors used Gaia DR2 sources with astrometric_excess_noise=0, only a few of which fall in the central region of the ONC, as shown in Figure 1. Hence, their sample does not fully reflect the distribution of stellar velocities over the region covered in this work. In Section 4.2, we identified deviation from normality at the tails of the distribution. To verify the multivariate normality more quantitatively, we employed the R package MVN (Korkmaz et al. 2014) based on the method of Henze & Zirkler (1990). As in Section 4.2, we tested the three cases: excluding (a) none of the ESC groups, (b) ESC group 1, and (c) ESC group 1 + 2. The first case exhibits statistically significant deviation from normality (p ∼ 0.01), while the latter two cases show consistency with a multivariate normal distribution (p > 0.05). The deviation from normality suggests that including the ESCs in ESC group 1 would overestimate the width of the PM distribution.

Understanding the nature of the kinematic outliers is important for assessing the applicability of our PM distribution model used in the previous section. We notice in Figure 7 that these stars are mostly concentrated at the cluster core (${r}_{c}\sim 0\buildrel{\,\prime}\over{.} 1;$ Hillenbrand & Hartmann 1998), with PM vectors heading outward. This is also the case for the stars in ESC group 2. Their positions and motions imply that the majority are likely higher-velocity stars escaping the ONC as a consequence of more frequent dynamical interactions between stars at the cluster core (Johnstone 1993; Baumgardt et al. 2002). Another possibility is that some of these are unresolved binaries centrally concentrated due to mass segregation, although in this case, their anisotropic radial PMs would not be easily explained. It is yet difficult to take a complete census of escaping stars due to measurement errors induced by the strong nebulosity and the lack of line-of-sight velocity measurements. With the relation $\langle {v}^{2}\rangle$ = ${\sigma }_{v}^{2}$ for the Maxwellian distribution and PM dispersions calculated in Section 4.3, Equation (1) would give an estimate for the mean-square escape speed of 2.8 ± 0.1 mas yr⁻¹, ∼10% lower than our previous estimate in Section 4.2. This suggests that the previous estimate for the escape velocity needs to be treated as an upper limit, and that applying the lower limit could reveal additional candidates. Nonetheless, if the higher-velocity stars are not escaping, then the deviation from normality seen in case (a) would be attributed to the rapid variation of velocity dispersions within the core of the cluster, as shown in Figure 8. We note that, even when including the ESCs, the PMs of all five subsamples grouped by distance in Table 5 are consistent with a multivariate normal distribution (p > 0.05) in both the Cartesian and radial-tangential coordinates. Our model is thus still valid in all of the radial bins.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 5.  PM Dispersions as a Function of Distance

Excluded	Radii	${\sigma }_{\mu ,{\alpha }^{* }}$	σμ,δ	σμ,r	σμ,t	${\sigma }_{\mu ,1{\rm{D}}}$	N
ESC Group
	arcmin	mas yr−1	mas yr−1	mas yr−1	mas yr−1	mas yr−1
	0.0−0.7	1.20 ± 0.10	1.40 ± 0.11	1.26 ± 0.10	1.41 ± 0.11	1.30 ± 0.07	91
	0.7−1.4	0.93 ± 0.06	1.18 ± 0.08	1.07 ± 0.07	1.10 ± 0.07	1.06 ± 0.05	153
None	1.4−2.1	0.80 ± 0.06	1.11 ± 0.08	0.96 ± 0.07	0.99 ± 0.07	0.97 ± 0.05	135
	2.1−2.8	0.80 ± 0.05	1.00 ± 0.06	0.99 ± 0.06	0.82 ± 0.05	0.91 ± 0.05	158
	2.8−3.5	0.81 ± 0.07	1.00 ± 0.08	0.83 ± 0.07	0.97 ± 0.08	0.91 ± 0.05	93
	0.0−0.7	0.99 ± 0.09	1.40 ± 0.11	1.13 ± 0.10	1.30 ± 0.11	1.21 ± 0.07	87
	0.7−1.4	0.90 ± 0.06	1.18 ± 0.08	1.02 ± 0.07	1.04 ± 0.07	1.04 ± 0.05	152
Group 1	1.4−2.1	0.80 ± 0.06	1.11 ± 0.07	0.96 ± 0.07	0.99 ± 0.07	0.97 ± 0.05	135
	2.1−2.8	0.80 ± 0.05	1.00 ± 0.06	0.99 ± 0.06	0.82 ± 0.05	0.91 ± 0.05	158
	2.8−3.5	0.75 ± 0.06	1.00 ± 0.08	0.84 ± 0.07	0.92 ± 0.08	0.88 ± 0.05	92
	0.0−0.7	0.98 ± 0.09	1.18 ± 0.10	1.04 ± 0.10	1.12 ± 0.10	1.08 ± 0.07	79
	0.7−1.4	0.90 ± 0.06	1.13 ± 0.07	0.99 ± 0.07	1.03 ± 0.07	1.02 ± 0.05	149
Group 1 + 2	1.4−2.1	0.80 ± 0.06	1.11 ± 0.07	0.96 ± 0.07	0.99 ± 0.07	0.97 ± 0.05	135
	2.1−2.8	0.80 ± 0.05	0.96 ± 0.06	0.96 ± 0.06	0.81 ± 0.05	0.88 ± 0.04	157
	2.8−3.5	0.75 ± 0.06	1.00 ± 0.08	0.84 ± 0.07	0.92 ± 0.08	0.88 ± 0.05	92

Download table as:  ASCIITypeset image

The ONC is well known to have a stellar velocity distribution that is elongated north–south (e.g., Jones & Walker 1988; Kuhn et al. 2019) along the axis of the Orion A cloud filament. The PM distribution of our sample also appears elongated north–south with an axis ratio ${\sigma }_{\mu ,{\alpha }^{* }}/{\sigma }_{\mu ,\delta }=0.74\pm 0.03$. This kinematic anisotropy agrees with that seen in the stellar distribution at a radius of ∼3', which is also elongated north–south with b/a ∼ 0.7 (see Figure 3 in Da Rio et al. 2014), parallel to the Orion A molecular cloud (Hillenbrand & Hartmann 1998). The PM dispersions may reflect the initial conditions of the protocluster cloud or the geometry of the present-day gravitational potential. On the other hand, the deviation from tangential to radial isotropy (${\sigma }_{\mu ,t}/{\sigma }_{\mu ,r}-1;$ Bellini et al. 2018) is 0.03 ± 0.04, which suggests that the cluster is consistent with being isotropic in the tangential–radial velocity space.

It has been argued in some previous studies that the ONC is likely to be supervirial (e.g., Jones & Walker 1988; Da Rio et al. 2014). The dynamical mass of the ONC inferred from the previous kinematic analysis done by Jones & Walker (1988) is nearly twice the total stellar and gas mass. Alternatively, virial equilibrium requires a one-dimensional mean velocity dispersion of σ_v ≃ 1.7 ± 0.3 km s⁻¹ given the volume density of the stellar and gas contents in the ONC (Da Rio et al. 2014), which is only ∼75% of the velocity dispersion of 2.34 ± 0.06 km s⁻¹ from Jones & Walker (1988), leading to the conclusion that the ONC is likely to be slightly supervirial, with a virial ratio of q ≃ 0.9 ± 0.3. This past result is, however, partially attributable to the previous estimate of the distance of the ONC adopted for the derivation of velocity dispersion in Jones & Walker (1988), ∼470 pc.

The estimate of the distance subsequently decreased in later studies to ∼400 pc (e.g., Jeffries 2007; Menten et al. 2007; Kounkel et al. 2017, 2018; Großschedl et al. 2018; Kuhn et al. 2019). At a distance of d = 414 ± 7 pc from Menten et al. (2007), Jones & Walker (1988) would have obtained a smaller value for their velocity dispersion, σ_v ≃ 2.1 ± 0.1 km s⁻¹, which would give a virial ratio as low as q ≃ 0.7 ± 0.3.

Figure 8 compares our measured velocity dispersions in Table 5 to the predicted velocity dispersions required for virial equilibrium based on the total mass profile from Da Rio et al. (2014). We note that we adopted 414 ± 7 pc from Menten et al. (2007) for the distance of the ONC, as Da Rio et al. (2014) did, for consistency purposes. Overall, our velocity dispersion profiles tend to decrease with radius following the prediction based on the observed mass profile, a power-law profile slightly steeper than a singular isothermal sphere (Da Rio et al. 2014). When we include the ESCs, our measured velocity dispersion appears to be more than 1σ larger than predicted at the very central region. However, when neglecting the escaping stars in ESC group 1 or 1 + 2, our measurements are in good agreement with the predicted values within 1σ uncertainty, suggesting that the bulk of the cluster is virialized. Even if the ESCs are included, the measured velocity dispersions are still well below the boundedness limit (${\sigma }_{\mathrm{bound}}=\sqrt{2}{\sigma }_{\mathrm{vir}}$). There is no indication of global expansion in the ONC from our analysis; the mean PM along the radial axis is small and consistent with zero within the σ uncertainty, as shown in Section 4.3. Kuhn et al. (2019) reported evidence of mild expansion in the ONC based on PM measurements of Gaia DR2 sources, finding the median outward velocity $\widetilde{{v}_{\mathrm{out}}}=0.42\pm 0.20$ km⁻¹ based on the uncertainty-weighted median with bootstrap resampling. In Figure 5 of their paper, the weighted kernel-density estimate (KDE) plot of v_out exhibits two peaks, one at ∼−0.4 km s⁻¹ and the other at ∼1.0 km s⁻¹, which implies a significant concentration of data points or weights at two different places. To explore this, we divided the Gaia DR2 sample of Kuhn et al. (2019), including 378 stars, into two subgroups in terms of weights: (a) 42 stars with small errors (i.e., large weights) _v,out < 0.15 km s⁻¹ and (b) 336 stars with larger errors _v,out > 0.15 km s⁻¹. The sum of the weights of group (a) reaches ∼64% of that of group (b). For groups (a) and (b), the median velocity is found to be $\widetilde{{v}_{\mathrm{out}}}=0.94\pm 0.41$ and 0.09 ± 0.19 km s⁻¹, respectively. We note that we converted PMs into v_out and calculated the median velocity in the same manner as described in Sections 3 and 4 of Kuhn et al. (2019), for consistency purposes. This result indicates that ∼10% of the whole sample (i.e., group (a)) produces the second peak in the KDE plot and biases the median velocity, ultimately leading to the conclusion that the ONC shows evidence for mild expansion. In fact, the stars in group (a) mostly fall into the magnitude range between phot_g_mean_mag ∼13 and 15, which corresponds to the transition point in the Gaia error terms from the detector and calibration-dominated regime to photon noise–limited regimes. In this interval, the astrometric uncertainties of Gaia DR2 are known to be the most underestimated (see Section 4.6.4 and Figure 24 in Arenou et al. 2018). We refer to the IAU GA presentation slides by Lindegren¹³ for more details. Between group (a) and our catalog, we found one star matched, whose PM errors in the Gaia DR2 need to be increased by a factor of ∼3 or greater for its outward velocities to be consistent within 1σ uncertainty (see also Section 4.1). We therefore conclude that the net outward velocity claimed by Kuhn et al. (2019) is likely a result of underestimated uncertainties in the Gaia DR2 data.

Our result adds to the growing evidence that the central region of the ONC is dynamically evolved. Studies of spatial morphology have revealed that the ONC has overall very little stellar substructure (Hillenbrand & Hartmann 1998; Da Rio et al. 2014). The core of the cluster exhibits a rounder and smoother stellar distribution than the outskirts, indicating that the core has likely experienced more dynamical timescales to lose initial substructures. The line-of-sight velocities are also smoothly distributed, as expected from an old dynamical age (Da Rio et al. 2017). Correcting for the local variation of the mean velocities, Da Rio et al. (2017) found that the dispersion of the line-of-sight velocities is measured as low as ∼1.7 km s⁻¹, which agrees with a virial state. Kuhn et al. (2019) also reached a similar conclusion using the PMs of 48 Gaia DR2 sources.

The dynamical age also has an implication for the origin of mass segregation. The ONC exhibits clear evidence of mass segregation, with the most massive stars preferentially located in the central region of the cluster (Hillenbrand & Hartmann 1998). The young age of the stellar population, as young as ∼2.5 Myr on average, is often cited as evidence for the primordial origin of the mass segregation. The evidence of virial equilibrium, however, implies that the central region of the ONC is dynamically old enough to have undergone several crossing times or more (Tan et al. 2006). With a large age spread of ∼1.3 Myr, the stellar population already appears several times older than the current crossing time based on the observed cluster mass (Da Rio et al. 2014). In addition, the dynamical timescale was likely much smaller at the early phase of the cluster due to higher stellar densities (Bastian et al. 2008; Allison et al. 2009). The mass segregation in the ONC thus need not be fully primordial.

Ultimately, our results favor the theoretical hypothesis that star cluster formation is a dynamically "slow" process; stars form slowly in supersonically turbulent gas over several crossing times, reaching a quasi-equilibrium (Tan et al. 2006; Krumholz & Tan 2007; Krumholz et al. 2012). In the "fast" scenarios, a star cluster forms via a rapid global collapse of a gas clump within approximately one crossing time (e.g., Elmegreen 2007; Krumholz et al. 2011), in which case the cluster would lack the features expected from dynamical evolution. The ONC continues to form stars with low efficiency and largely in virial equilibrium; thus, it provides observational evidence for the slow formation scenario and does not support the competitive accretion model.

Our PM catalog includes previously known fast-moving sources around the BN/KL complex, namely, BN, source x, and source n, although in our analysis, source n exhibits a rather small PM in the rest frame of the ONC. These optically invisible objects are all detected in epoch 2010 and 2014 NIRC2 He i b (∼K-band) images.¹⁴ Before the recent discovery of source x, two possible scenarios were proposed to explain their origin:

(a) Ejected from Trapezium θ¹ Ori C ∼4000 yr ago (Plambeck et al. 1995; Tan 2004), BN passed near source I, triggering an explosive outflow; or (b) BN, source I, and at least one other star once comprised a multiple system, and the latter two objects merged into a tight binary system, which resulted in an explosive outflow ∼500 yr ago (Rodríguez et al. 2017).

Numerous pieces of evidence argue against the first scenario, e.g., the large separation between θ¹ Ori C and BN at the time of ejection (≳10'') and their inconsistent momenta and ages (see Goddi et al. 2011 and references therein). In the meantime, source x was recently found as a promising candidate of the formerly missing puzzle for the second scenario. Luhman et al. (2017) demonstrated that the estimated location of source x ∼500 yr ago agrees with the position for BN and source I at that time, determined by radio PMs measured by Gómez et al. (2008), Goddi et al. (2011), and Rodríguez et al. (2017), which suggests that the three sources were likely ejected in the same event.

We have independently calculated when BN and source x experienced their closest approach (t_min) using our PM measurements and the method described in Appendix B. We obtained t_min = 1535 ± 29 yr, which is consistent within 2σ with that for BN and source I based on their radio PMs, t_min = 1475 ± 6 yr (Rodríguez et al. 2017). Figure 9 shows the PM vectors and positions of BN and sources I, n, and x and the 1σ range of allowed paths back to 1535, where only the data for source I have been adopted from Rodríguez et al. (2017). Our observation supports the possibility that BN and sources I and x were ejected from around the same location ∼500 yr ago. Elements of this hypothesis were once challenged by theoretical simulations; Farias & Tan (2018) demonstrated using a large suite of N-body simulations that the mass for source I required by their simulations is as large as at least 14 M_⊙, almost twice as large as the previously estimated ∼7 M_⊙ based on the kinematics of the circumstellar material (Matthews et al. 2010; Hirota et al. 2014; Plambeck & Wright 2016). Recent ALMA observations with higher resolution and sensitivity than the previous ones, however, estimate the mass of source I as $15\pm 2{M}_{\odot }$ (Ginsburg et al. 2018), by which the dynamical decay scenario still remains viable.

Zoom In Zoom Out Reset image size