The 3D Kinematics of the Orion Nebula Cluster: NIRSPEC-AO Radial Velocities of the Core Population

All objects in our curated sample were observed using NIRSPEC on Keck II, in conjunction with the laser guide star (LGS) AO system (van Dam et al. 2006; Wizinowich et al. 2006). We used the instrument in its high spectral resolution AO mode with a slit width × slit length of 0041 × 226. Observations were carried out in the K or N7 band to capitalize on the strong CO bands within that regime (∼2.3 μm, orders 32 and 33). We additionally used a cross-disperser angle of 3565 and an echelle angle of 63°. The resolution in this setup is R ≈ 25,000, as determined by the width of unresolved OH sky lines, and the wavelengths covered are approximately 2.044–2.075 (order 37), 2.100–2.133 (order 36), 2.160–2.193 (order 35), 2.224–2.256 (order 34), 2.291–2.325 (order 33), and 2.362–2.382 (order 32) μm, with some portions of the bands beyond the edges of the detector. For this work, all analysis was done using orders 32 and 33, the orders containing the CO bandheads in addition to numerous telluric features, and all NIRSPAO data presented come from these orders. We note that NIRSPEC underwent an upgrade during 2018, and our setup changed slightly postupgrade. The most significant change is a larger wavelength coverage for each order due to a larger Hawaii 2RG detector (2048 × 2048 postupgrade versus 1024 × 1024 preupgrade; Martin et al. 2018) and a higher resolution (R ≈ 35,000; Hsu et al. 2021b).

Typical observations consisted of four spectra taken in an ABBA dither pattern along the length of the slit. In a few cases, more or less than four spectra were taken. Either before or after each target was observed, an A0V calibrator star at similar airmass was observed for a telluric reference. Table 1 gives the log of our spectroscopic observations, listing the targets observed, the date of observation, the number of spectra, and the integration time for each spectrum. Etalon lamp and flat-field frames were also taken each night for use in data reduction (Section 3).

Our data were forward-modeled using the Spectral Modeling Analysis and RV Tool (SMART; ¹⁰ Hsu et al. 2021a). Details for the fitting routine using SMART are given in Hsu et al. (2021b). Here we briefly outline our methods.

Reduced NIRSPAO data were modeled using an iterative approach. The first step was obtaining an absolute wavelength solution to orders 32 and 33. For our initial wavelength solution, we use the global wavelength solution provided by the NSDRP, which is a quadratic polynomial ¹¹ of the form

$$\begin{eqnarray}&&\lambda (p,M)={r}_{0}+{r}_{1}p+{r}_{2}{p}^{2}+\displaystyle \frac{{r}_{3}}{M}+\displaystyle \frac{{r}_{4}p}{M}+\displaystyle \frac{{r}_{5}{p}^{2}}{M},\end{eqnarray} \tag{ 1 }$$

where λ(p, M) is the wavelength falling on column pixel p of order M and coefficients r_n . This initial wavelength solution is obtained from fitting to the etalon lamps, which provides uniform offsets in frequency space. However, as mentioned previously, the absolute wavelength solution, or starting position, is unknown for the etalon spectrum.

To obtain a more precise wavelength calibration, we performed a cross-correlation between our telluric spectrum from the A-star calibrator to the high-resolution telluric spectrum from Moehler et al. (2014). First, we modeled and removed the continuum of our A-star calibrator using a quadratic polynomial, essentially leaving just the flat imprinted telluric absorption spectrum. Then, we scaled the flux of the high-resolution telluric model to the A-star flux. Next, the telluric spectrum from the A0V star was cross-correlated with the telluric model using a window of 100 pixels and a step size of 20 pixels, calculating the best-fit cross-correlation shift for each window along the entire spectrum. Then, a fourth-order polynomial—i.e., λ(p) = a_i + b_i p + c_i p² +d_i p³ + e_i p⁴, where i is the iteration number—was fit to the best wavelength shifts for all windows used in that iteration. The initial wavelength solution was given using the second-order polynomial provided by the NSDRP, with the additional coefficients set to zero (i.e., d₀, e₀). The best-fit coefficients for each pass (e.g., ${\delta }_{{a}_{0}}$) were added to the previous solution, e.g., ${a}_{1}={a}_{0}+{\delta }_{{a}_{1}}$, ${b}_{1}={b}_{0}+{\delta }_{{b}_{1}}$. This loop was repeated until the wavelength solution converged to the smallest residuals between the telluric spectrum and telluric model. One thing to note is that different iterations used a different pixel window and step size, as finer granularity is required as the wavelength solution gets closer to the optimal fit. For instance, the second pass used a step size of 10 pixels and a window size of 150 pixels.

The aforementioned fitting procedure was done for each frame independently, resulting in an absolute wavelength solution for each frame. We cross-correlated a single A-star frame to A-star calibration data taken over 14 nights (both A and B nods) and found an rms between frames of 0.004 Å (0.058 km s⁻¹). Although this systematic uncertainty is typically much smaller than our measurement uncertainty, we add this systematic in quadrature with our measurement uncertainties per frame.

Next, we modeled the spectrum using a forward-modeling approach based on the method provided by Blake et al. (2010) and also Butler et al. (1996), Blake et al. (2007, 2008, 2010), and Burgasser et al. (2016), utilizing a Markov Chain Monte Carlo (MCMC) method built on emcee (Foreman-Mackey et al.2013) to sample the parameter space.

The flux from the source can be modeled using the following equation:

$$\begin{eqnarray}\begin{array}{rcl}{F}_{M}[p] & = & C[p(\lambda )]\times \left[\left(\left(M\left[{p}^{\ast }\left(\lambda \left[1+\displaystyle \frac{{V}_{r}}{c}\right]\right),{T}_{{\rm{eff}}},{\rm{log}}g,[{\rm{M}}/{\rm{H}}]\right])\right.\right.\right.\\ & & \ast \left.\left.\left.{\kappa }_{R}(v\sin i)+{C}_{{\rm{veil}}}\right)\times T[{p}^{\ast }(\lambda ),{\rm{AM}},{\rm{PWV}}\right]\right]\\ & & \ast {\kappa }_{G}({\rm{\Delta }}{\nu }_{{\rm{inst}}})+{C}_{{\rm{flux}}},\end{array}\end{eqnarray} \tag{ 2 }$$

where an asterisk indicates convolution, C[p(λ)] is a quadratic polynomial that is used for continuum correction (meant to correct and scale for variations induced by the instrumental profile on the observed and absolute flux of the model), and M is the photospheric model of the source. We fixed $\mathrm{log}g=4$, as this is approximately the expected gravity for low-mass stars at the age of the ONC that are still contracting onto the main sequence and consistent with other studies (e.g., Kounkel et al. 2018). Additionally, we show later that RV variations with $\mathrm{log}g$ tend to be small. We also fixed [M/H] = 0 based on the average metallicity of the ONC (e.g., D'Orazi et al. 2009). Here ${\kappa }_{R}(v\sin i)$ is the rotational broadening kernel (using the methods of Gray 1992 and a limb-darkening coefficient of 0.6); V_r and c are the heliocentric RV and speed of light, respectively; T is the telluric spectrum from Moehler et al. (2014); and AM and PWV are the airmass and precipitable water vapor of the telluric spectrum, respectively. Veiling is parameterized by C_veil, which is an additive graybody flux to the stellar model flux (continuum), r_λ = F_λ,cont/F_λ,* (or C_veil/F_λ,*), to represent potential veiling along the line of sight, which will weaken the depths of the photospheric absorption lines (e.g., Muzerolle et al. 2003a, 2003b; Fischer et al. 2011). Extinction effects, which reduce the intensity of the emission lines, are multiplicative in nature and therefore get folded into the fit for C[p(λ)]. We note that due to the small range of wavelengths used, effects due to extinction should be minimal and would mostly impact stellar parameters, rather than RVs.

Here p* = p(λ) + C_λ , where p(λ) is a fourth-order polynomial mapping of pixel to wavelength based on the telluric spectrum from the absolute wavelength calibration to the A star, and C_λ is a small constant offset/correction to the zeroth-order term. We keep constant all coefficients in the fourth-order polynomial that were derived from the A-star calibrator, save for the zeroth-order term, which we fit for a small constant offset (nuisance parameter) to account for small differences in the absolute wavelength calibration versus the observed data. This is necessary, as observations not taken along the exact same pixels will shift by a small amount due to the spatial curvature of the flux along the detector. Here κ_G (Δν_inst) is the spectrograph line-spread function (LSF), modeled as a normalized Gaussian. We include an additive flux offset, C_flux, as an additional nuisance parameter to account for small differences in the absolute flux calibration. We fit for separate C_λ and C_flux parameters for each order. For stellar models, we chose the PHOENIX-ACES-AGSS-COND-2011 stellar models (Husser et al. 2013), which have been used previously for modeling ONC young stellar objects (YSOs) observed in the NIR with SDSS/APOGEE (Kounkel et al. 2018)

The log-likelihood function, assuming normally distributed parameters and noise, is

$$\begin{eqnarray}&&\begin{array}{c}{\rm{l}}{\rm{n}}{ \mathcal L }=-\displaystyle \frac{1}{2}\sum _{p}\left[{\left(\displaystyle \frac{D[p]-{F}_{M}[p]}{\sigma [p]\times {C}_{{\rm{n}}{\rm{o}}{\rm{i}}{\rm{s}}{\rm{e}}}}\right)}^{2}\right.\\ \,+\,\left.{\rm{l}}{\rm{n}}[2\pi {(\sigma [p]\times {C}_{{\rm{n}}{\rm{o}}{\rm{i}}{\rm{s}}{\rm{e}}})}^{2}]\right],\end{array}\end{eqnarray} \tag{ 3 }$$

where D[p] is the data, σ[p] is the noise or flux uncertainty, and C_noise is a scaling parameter for the noise to account for systematic errors between the observations and the models.

We simultaneously fit for all of the above parameters, using the affine-invariant ensemble sampler emcee, with the kernel-density estimator described in Farr et al. (2014). Uniform priors were used across the parameter ranges shown in Table 2. We did an initial fit using 100 walkers and 400 steps, discarding the first 300 steps. Typical convergence occurred after the initial 200 steps. These fits were then masked for bad pixels outside of three standard deviations of the median difference between the model and the data, effectively removing bad pixels and cosmic rays that were not removed by the fixpix utility. Next, another fit was performed on the masked data with 100 walkers, 300 steps, and a burn-in of 200 steps. Walkers were initialized within 10% of the best-fit parameters from the initial fit. Convergence typically occurred within 100 steps. Heliocentric RVs were corrected for barycentric motion using the astropy function radial_velocity_correction.

Table 2. Forward-modeled Parameter Ranges

Description	Symbol	Bounds
Stellar effective temp.	Teff	(2300, 7000) K
Rotational velocity	$v\sin i$	(1, 100) km s−1
RV	Vr	(−100, 100) km s−1
Airmass	AM	(1, 3)
Precip. water vapor	PWV	(0.5, 30) mm
LSF	Δνinst	(1, 100) km s−1
Flux offset param.	Cflux	(10−15, 1015) counts s−1
Noise factor	Cnoise	(1, 50)
Wave. offset param.	Cλ	(−10, 10) Å

Download table as:  ASCIITypeset image

An example fit is shown in Figure 2, with the observed—telluric-corrected—spectrum shown with the gray line, the best-fit stellar model shown with the red line (model parameters indicated in red text), and the best-fit stellar model convolved with the best-fit telluric model shown with the purple line. The gray spectrum and purple model are compared in our MCMC routine. The bottom plots show the residuals between the observed spectrum and best-fit model (black line) and the total uncertainty (noise computed from the NSDRP with the scaling factor in Equation (3)). The largest contributor to the residuals is fringing, which is an ongoing project to model and will be addressed in a future study. However, previous studies of modeling NIRSPEC fringing have shown that it does not significantly effect RV determinations (Blake et al. 2010). Figure 3 shows the corner plot for the MCMC run for Figure 2 (last 100 steps × 100 walkers, 10,000 data points). The parameters listed on the x- and y-axes are the same as those from Equations (2) and (3). This plot indicates that there is little to no correlation between the majority of parameters, with the exception of $v\sin i$ and veiling, which show a slight correlation. As can be seen, some of these distributions are non-Gaussian, which is why we report median values with 16th and 84th percentiles as uncertainties. The results of our fits are listed in Table 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Many of our sources converged to relatively high veiling parameters (i.e., r_λ > 0.2), and we note that temperature should be highly degenerate with veiling ratio due to the strong dependence of the fits on the CO bands, which could be weakened either by higher effective temperatures or higher veiling ratios. However, we do not have sufficient data to put constraints on the veiling parameter for each source, as that would required a 3D extinction map of the ONC (e.g., Schlafly et al. 2015) and precise distances to individual sources. However, even with the Gaia satellite (Gaia Collaboration et al. 2016), the embedded nature of the ONC makes parallax measurements extremely difficult (Kim et al. 2019). Therefore, this study is focused on the kinematics of the ONC; however, a future study will investigate the T_eff (and mass) dependence of kinematics in the ONC core.

Our RV measurements are relatively robust to changes in stellar parameters, since they are strongly anchored to the CO bandheads and the absolute calibration of the telluric spectrum. To illustrate this, we fit [HC 2000] 244 (the same source shown in Figure 2) using the same procedure outlined above and holding the $\mathrm{log}g$ constant from zero to 6 in steps of 0.5 (the resolution of the model grid). Figure 4 (top) shows the RV variation due to different $\mathrm{log}g$ values. For $\mathrm{log}g$ values that differ by more than ±1 dex from the nominal value, RV variations are more than 1%; however, those are extremely large variations in $\mathrm{log}g$ that are inconsistent with the youth of these targets. For RV values within ±0.5 dex, the variation is less than 0.5%. To be conservative, we chose to add this systematic uncertainty to our combined measurement uncertainty across all frames for each single object.

Zoom In Zoom Out Reset image size

We performed a similar test to the one above, this time holding temperature constant. Figure 4 (bottom) shows the RV variation due to different temperatures. Within a few hundred kelvins, the RV variation is less than 3%. However, within ±100 K, this value is less than 0.5%, which is within the systematic uncertainty found in the $\mathrm{log}g$ variation. Therefore, no additional systematic is required, since the ±0.5% systematic accounts for both the $\mathrm{log}g$ and T_eff variations. In summary, our reported uncertainties are the summed quadrature of our measurement uncertainty from the MCMC fits, the 0.058 km s⁻¹ systematic uncertainty between calibration frames, and the 0.5% variation found from differing $\mathrm{log}g$ and T_eff.

Lastly, we compare our derived RVs to those from APOGEE (using the results of K18) and Tobin et al. (2009), using the reanalyzed RVs from Kounkel et al. (2016, hereafter K16). Figure 5 shows the results of our RV comparison. In total, there were 11 matches between our NIRSPEC targets and the T09/K16 sample and seven matches to APOGEE (using a 05 crossmatch radius). Our results are consistent with the K18 APOGEE results, with only two exceptions differentiating by more than 1σ of their combined uncertainty, possibly the result of spectroscopic binaries: 2MASS J05352104−0523490 and 2MASS J05351498−0521598. In comparison to optical RVs from T09/K16, our measured RVs are generally smaller than those measured from optical data, by ∼1.8 km s⁻¹, on average. We also compared K16 to K18 RVs, again using a 05 crossmatch radius, finding 586 sources in common (Figure 5, green symbols). The distribution of the uncertainty-weighted difference between these two measurements, $({\mathrm{RV}}_{{\rm{K}}18}-{\mathrm{RV}}_{{\rm{K}}16})/\sqrt{{\sigma }_{{\mathrm{RV}}_{{\rm{K}}18}}^{2}+{\sigma }_{{\mathrm{RV}}_{{\rm{K}}16}}^{2}}$, has a mean value of μ = 0.58 km s¹ and a standard deviation of σ = 1.83 km s¹. This is consistent with the K16 values being, on average, smaller than the K18 APOGEE RVs. The width of the distribution also indicates that one or both of the uncertainties are underestimated. In general, RVs derived from NIR data versus optical are likely to suffer from fewer systematics in this highly embedded region.

Zoom In Zoom Out Reset image size

Our measurements provide velocities along the line of sight; however, to measure 3D velocities, we require tangential motions. A number of studies have measured the PMs of sources within the ONC (e.g., Parenago 1954; Jones & Walker 1988; van Altena et al. 1988; Gómez et al. 2005; Dzib et al. 2017; Kim et al. 2019; Kuhn et al. 2019; Platais et al. 2020). The two most recent catalogs produced by Kim et al. (2019) and Platais et al. (2020) both use imaging data from the Hubble Space Telescope (HST). We chose to adopt the values from Kim et al. (2019), as their combination of HST imaging with data from the Keck Near Infrared Camera 2 (NIRC2; PI: K. Matthews) provides a baseline of ∼20 yr. Additionally, the uncertainties published by Platais et al. (2020) tend to be very small and are possibly underestimated due to the fact that they were unable to determine systematic uncertainties. This likely explains their much smaller errors versus Kim et al. (2019), even though a shorter time baseline of ∼11 yr was used.

The Kim et al. (2019) catalog includes PM measurements for 701 sources with typical uncertainties ≲1 mas yr⁻¹. However, to convert PMs into tangential velocities, we require a distance, or distance distribution, to the ONC. A number of studies have investigated the distance to the ONC using very long baseline interferometry (VLBI; Menten et al. 2007; Sandstrom et al. 2007; Kounkel et al. 2017) and, more recently, the Gaia Data Release 2 (DR2; Gaia Collaboration et al. 2018; Großschedl et al. 2018; Kounkel et al. 2018; Kuhn et al. 2019). The majority of these measurements are consistent with an average distance to the ONC of ∼390 pc. For our study, we adopted the VLBI trigonometric distance of 388 ± 5 pc from Kounkel et al. (2017), which is consistent with Kounkel et al. (2018; 389 ± 3 pc, Gaia DR2), Kuhn et al. (2019; ${403}_{-6}^{+7}$ pc, Gaia DR2), Großschedl et al. (2018; 397 ± 16 pc, Gaia DR2), and Sandstrom et al. (2007; ${389}_{-21}^{+24}$ pc, VLBI).

Using the above distance estimate, we converted PMs to tangential velocities, combining errors in PMs and the distance to the ONC using standard error propagation. The combined sample with all three components of motion totaled 135 sources. Figure 8 shows a map of the ONC with vectors displaying the PM of sources in our sample and colors representing the measured RVs. All three components of motion for our subsample are shown in Figure 9. These velocities were used to investigate the 3D kinematics of the ONC core region. There are a number of kinematic outliers, both in tangential velocity space (as discussed in Kim et al. 2019) and in RV space, that will be discussed in Section 5.3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To determine velocity dispersion, we utilized a similar Bayesian framework to Kim et al. (2019), where each ith kinematic measurement is parameterized as ${v}_{{\left(\alpha ,\delta ,r\right)}_{i}}\pm {\epsilon }_{{\left(\alpha ,\delta ,r\right)}_{i}}$. We assume each kinematic measurement is drawn from a multivariate Gaussian distribution with mean values of ${\bar{v}}_{\alpha }$, ${\bar{v}}_{\delta }$, and ${\bar{v}}_{r}$ and standard deviations of $\sqrt{{\sigma }_{{v}_{\left(\alpha ,\delta ,r\right)}}^{2}+{\epsilon }_{{\left(\alpha ,\delta ,r\right)}_{i}}^{2}}$, where ${\sigma }_{{v}_{\alpha }}$, ${\sigma }_{{v}_{\delta }}$, and ${\sigma }_{{v}_{r}}$ are the intrinsic velocity dispersions (IVDs). The log-likelihood for the 3D kinematics of the ith object is given by

$$\begin{eqnarray}&&\mathrm{ln}{L}_{i}=-\displaystyle \frac{1}{2}\left[\mathrm{ln}(| {{\boldsymbol{\Sigma }}}_{i}| )+{\left({{\boldsymbol{v}}}_{i}-{\boldsymbol{\mu }}\right)}^{{\mathsf{T}}}{{\boldsymbol{\Sigma }}}_{i}^{-1}({{\boldsymbol{v}}}_{i}-{\boldsymbol{\mu }})+3\mathrm{ln}(2\pi )\right],\end{eqnarray} \tag{ 4 }$$

where the input kinematic measurement vector for the ith object v _i is defined as

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{i}=\left[\begin{array}{c}{v}_{{\alpha }_{i}}\\ {v}_{{\delta }_{i}}\\ {v}_{{r}_{i}}\end{array}\right],\end{eqnarray} \tag{ 5 }$$

while the vectors for the mean velocities and dispersions are defined as

$$\begin{eqnarray}&&{\boldsymbol{\mu }}=\left[\begin{array}{c}{\bar{v}}_{\alpha }\\ {\bar{v}}_{\delta }\\ {\bar{v}}_{r}\end{array}\right],\qquad {\boldsymbol{\sigma }}=\left[\begin{array}{c}{\sigma }_{{v}_{\alpha }}\\ {\sigma }_{{v}_{\delta }}\\ {\sigma }_{{v}_{r}}\end{array}\right],\end{eqnarray} \tag{ 6 }$$

and the covariance matrix Σ_i for the ith object is defined as

$$\begin{eqnarray}&&\begin{array}{l}{{\boldsymbol{\Sigma }}}_{i}\\ \,=\left[\begin{array}{ccc}{\sigma }_{{v}_{\alpha }}^{2}+{\epsilon }_{{\alpha }_{i}}^{2} & {\rho }_{1}({\sigma }_{{v}_{\alpha }}{\sigma }_{{v}_{\delta }}+{\epsilon }_{{\alpha }_{i}}{\epsilon }_{{\delta }_{i}}) & {\rho }_{2}({\sigma }_{{v}_{\alpha }}{\sigma }_{{v}_{r}}+{\epsilon }_{{\alpha }_{i}}{\epsilon }_{{r}_{i}})\\ {\rho }_{1}({\sigma }_{{v}_{\alpha }}{\sigma }_{{v}_{\delta }}+{\epsilon }_{{\alpha }_{i}}{\epsilon }_{{\delta }_{i}}) & {\sigma }_{{v}_{\delta }}^{2}+{\epsilon }_{{\delta }_{i}}^{2} & {\rho }_{3}({\sigma }_{{v}_{\delta }}{\sigma }_{{v}_{r}}+{\epsilon }_{{\delta }_{i}}{\epsilon }_{{r}_{i}})\\ {\rho }_{2}({\sigma }_{{v}_{\alpha }}{\sigma }_{{v}_{r}}+{\epsilon }_{{\alpha }_{i}}{\epsilon }_{{r}_{i}}) & {\rho }_{3}({\sigma }_{{v}_{\delta }}{\sigma }_{{v}_{r}}+{\epsilon }_{{\delta }_{i}}{\epsilon }_{{r}_{i}}) & {\sigma }_{{v}_{r}}^{2}+{\epsilon }_{{r}_{i}}^{2}\end{array}\right].\end{array}\end{eqnarray} \tag{ 7 }$$

The correlation coefficients are given by ρ_1,2,3 and should be zero if the covariance is completely uncorrelated. Using Bayes's theorem, for a set of N measurements ${\boldsymbol{D}}\equiv {\{{{\boldsymbol{v}}}_{i},{{\boldsymbol{\epsilon }}}_{i}\}}_{i=1}^{N}$, the log of the posterior probability is given as

$$\begin{eqnarray}&&\mathrm{ln}P({\boldsymbol{\mu }},{\boldsymbol{\sigma }}| {\boldsymbol{D}})\propto \sum _{i}^{N}\left[\mathrm{ln}{L}_{i}({{\boldsymbol{v}}}_{i},{{\boldsymbol{\epsilon }}}_{i}| {\boldsymbol{\mu }},{\boldsymbol{\sigma }})+\mathrm{ln}p({\boldsymbol{\mu }},{\boldsymbol{\sigma }})\right],\end{eqnarray} \tag{ 8 }$$

where p( μ , σ ) ≡ p( μ )p( σ ) is the prior on the means and dispersions. Similar to Kim et al. (2019), we adopted a flat prior on the mean vector, μ , and a Jeffreys prior on the dispersion vector, p( σ ) ∝ σ ⁻¹. To determine best-fit parameters for ${\bar{v}}_{\alpha }$, ${\bar{v}}_{\delta }$, ${\bar{v}}_{r}$, ${\sigma }_{{v}_{\alpha }}$, ${\sigma }_{{v}_{\delta }}$, ${\sigma }_{{v}_{r}}$, ρ₁, ρ₂, and ρ₃, we utilized the MCMC affine-invariant ensemble sampler emcee (Foreman-Mackey et al.2013).

Previous studies have identified subpopulations within the ONC with kinematics consistent with high-velocity escaping stars (e.g., Kim et al. 2019; Kuhn et al. 2019; Platais et al. 2020). Kim et al. (2019) and Kuhn et al. (2019) used outliers from an expected normal distribution to determine kinematic outliers. Kim et al. (2019) also used the ONC's 2D mean-square escape velocity (≈3.1 mas yr⁻¹, or 6.1 km s⁻¹ at a distance of 414 pc) to identify potential escaping sources. It should be noted that Kim et al. (2019) used a distance of 414 ± 7 pc from Menten et al. (2007), whereas our adopted value of 388 ± 5 pc would give a mean-square escape velocity of 5.7 km s⁻¹. In practice, this should not alter the results, since outliers are computed from the bulk properties of the cluster, and our smaller distance would equally impact all tangential velocities the same.

To determine high-probability escaping or evaporating sources, we identified sources whose velocities deviate from the Gaussian velocity distribution model (see Section 5.2). We used a methodology similar to Kim et al. (2019) and Kuhn et al. (2019), plotting sources on Q–Q plots, where data quantiles (Q_data) are plotted against theoretical quantiles (Q_theo) corresponding to a Gaussian distribution. These two quantiles are defined as

$$\begin{eqnarray}&&{Q}_{\mathrm{data},i}=\displaystyle \frac{{v}_{{\left(\alpha ,\delta ,r\right)}_{i}}-{\bar{v}}_{(\alpha ,\delta ,r)}}{\sqrt{{\sigma }_{\left(\alpha ,\delta ,r\right)}^{2}+{\epsilon }_{{\left(\alpha ,\delta ,r\right)}_{i}}^{2}}},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{Q}_{\mathrm{theo},i}=\sqrt{2}\ {\mathrm{erf}}^{-1}\left(\displaystyle \frac{2({r}_{i}-0.5)}{n}-1\right),\end{eqnarray} \tag{ 10 }$$

where erf⁻¹ is the inverse of the error function, r_i is the rank of the ith measurement, and n is the number of measurements.

The means and IVDs are computed using the method outlined in Section 5.2. Figure 10 shows the quantiles for all sources (top panels), quantiles after removing outliers in RV space (middle panels), and quantiles after removing large outliers in RV and tangential kinematic space (bottom panels). In each panel, the gray band indicates the 95% confidence interval from bootstrapping. Kim et al. (2019) identified two separate groups of kinematic outliers, sources with highly significant deviations within the Q–Q plot (Q_data,X ≥ ±3; escape group 1) and sources with low-significance escape velocities defined from the angular escape speed (sources with μ_tot > 3.1 mas yr⁻¹; escape group 2). We additionally identify sources with RVs that significantly deviate from the sample (Q_data,Z ≥ ±3) and label this escape group 3 (see Figures 9 and 10). Kim et al. (2019) and Kuhn et al. (2019) estimated the escape speed for the ONC to be ≈6.1 km s⁻¹ using the virial theorem, and all of the outliers in escape groups 1 and 3 have velocities in excess of this speed. However, it should be noted that velocity outliers in RV space could be potential binaries rather than escaping/evaporating sources.

Zoom In Zoom Out Reset image size

Of the escaping stars identified in Kim et al. (2019), one of the sources in escape group 1 and three of the sources in escape group 2 were observed with NIRSPAO, and none were in our reanalyzed APOGEE subsample (identified in Table 3). This is not surprising, as all of the sources in escape groups 1 and 2 are within 25 of the ONC CoM (16 of 18 sources within 1'). The majority of sources in escape groups 1 and 2 tend to be fainter (14 of 18 sources with F139 mag ≳12) or have a nearby source within 2''. All of the sources in escape groups 1 and 2 have RVs that are consistent with the bulk RV distribution of the ONC. Conversely, none of the sources in escape group 3 (line-of-sight velocity outliers) are within escape groups 1 and 2, and both escape group 3 sources are APOGEE sources. This indicates that the high-velocity components for these sources are primarily in either the tangential or line-of-sight direction, but not both.

Table 3. NIRSPEC Forward-modeling Results

(HC 2000)	ID a	αJ2000	δJ2000	$\overline{{RV}}$ b	${\mu }_{\alpha \cos \delta }$	μδ	Teff c	$v\sin i$	Veiling d	Note e
		(deg)	(deg)	(km s−1)	(mas yr−1)	(mas yr−1)	(K)	(km s−1)	$\left(\displaystyle \frac{{F}_{{\rm{O}}33,\mathrm{cont}}}{{F}_{{\rm{O}}33,* }}\right)$
135	⋯	83.80950000	−5.40686111	28.56 ± 0.91	⋯	⋯	3610 ± 300	49.33 ± 2.26	0.24
188	559	83.83175000	−5.39925000	29.74 ± 0.27	0.41 ± 0.52	−0.83 ± 0.18	3475 ± 73	6.90 ± 1.72	0.01
202	615	83.81666667	−5.39805556	25.69 ± 0.75	2.80 ± 0.06	−1.56 ± 0.08	3393 ± 204	42.27 ± 2.95	0.75	E1
204	⋯	83.84087500	−5.39830556	24.73 ± 0.55	⋯	⋯	3712 ± 50	48.16 ± 2.56	0.28
215	⋯	83.80116667	−5.39669444	27.91 ± 0.27	⋯	⋯	3641 ± 33	9.03 ± 1.41	0.20
217	⋯	83.83770833	−5.39694444	29.73 ± 0.35	⋯	⋯	3614 ± 33	47.35 ± 0.89	0.23	V
219	87	83.81316667	−5.39630556	24.08 ± 0.48	1.87 ± 0.06	−0.65 ± 0.02	3486 ± 57	20.78 ± 1.40	0.43
220	23	83.81175000	−5.39625000	33.47 ± 0.23	−1.66 ± 0.40	0.08 ± 0.13	3484 ± 19	9.04 ± 0.30	0.02
223	164	83.81433333	−5.39597222	27.45 ± 0.23	1.10 ± 0.08	−0.68 ± 0.23	3359 ± 16	9.74 ± 0.69	0.32
224	⋯	83.79475000	−5.39575000	25.94 ± 0.15	⋯	⋯	3859 ± 21	24.94 ± 0.58	0.00
228	⋯	83.80250000	−5.39561111	27.46 ± 0.27	⋯	⋯	3782 ± 30	40.37 ± 0.78	0.46
229	536	83.83904167	−5.39594444	28.30 ± 0.21	0.12 ± 0.14	0.67 ± 0.27	3613 ± 79	12.19 ± 1.21	0.15
240	⋯	83.80604167	−5.39455556	27.11 ± 0.44	⋯	⋯	3640 ± 65	14.15 ± 1.28	3.89
244	180	83.82108333	−5.39438889	28.18 ± 0.11	−0.02 ± 0.15	0.83 ± 0.40	3436 ± 17	21.29 ± 0.83	0.04
245	⋯	83.81229167	−5.39425000	31.05 ± 0.21	⋯	⋯	3869 ± 18	17.03 ± 0.35	0.00
248	200	83.81566667	−5.39400000	26.56 ± 0.78	1.81 ± 0.43	−2.83 ± 0.20	3438 ± 48	39.39 ± 3.82	0.35	E2
250	197	83.81625000	−5.39388889	29.93 ± 0.46	1.70 ± 0.40	−3.85 ± 1.18	2975 ± 91	15.90 ± 1.86	0.39	E2
253	⋯	83.82587500	−5.39330556	28.67 ± 0.36	⋯	⋯	4000 ± 117	6.21 ± 1.48	2.66
258	521	83.81000000	−5.39269444	28.02 ± 0.19	−0.11 ± 0.23	−1.65 ± 0.03	3568 ± 41	7.71 ± 1.03	0.82
259	⋯	83.82112500	−5.39277778	27.53 ± 0.08	⋯	⋯	3501 ± 44	32.49 ± 1.95	0.54
261	206	83.81408333	−5.39261111	26.39 ± 0.20	−0.72 ± 0.24	0.87 ± 0.12	3374 ± 14	1.49 ± 0.31	0.00
275	65	83.81220833	−5.39141667	30.88 ± 0.01	−1.18 ± 0.14	1.19 ± 0.06	3860 ± 17	17.59 ± 0.30	0.00
277A	530	83.83525000	−5.39158333	25.88 ± 0.36	−0.81 ± 0.12	0.14 ± 0.06	3444 ± 51	12.90 ± 2.65	0.36	B
282	44	83.82866667	−5.39136111	26.96 ± 0.28	−0.85 ± 0.24	1.35 ± 0.03	3394 ± 22	5.59 ± 1.05	1.06
288	71	83.82966667	−5.39086111	26.61 ± 0.17	−0.42 ± 0.04	−0.52 ± 0.41	3666 ± 54	18.20 ± 0.53	0.25
291	211	83.81600000	−5.39044444	29.06 ± 0.22	1.10 ± 0.45	−1.63 ± 0.07	3181 ± 33	15.83 ± 0.72	0.42	B
295	⋯	83.82320833	−5.39025000	21.85 ± 0.93	⋯	⋯	3662 ± 305	22.69 ± 7.18	0.76
302	221	83.81133333	−5.38969444	29.24 ± 0.11	−0.24 ± 0.40	0.71 ± 0.13	3541 ± 21	13.73 ± 1.12	0.48
306A	⋯	83.81600000	−5.38958333	29.12 ± 0.50	⋯	⋯	4201 ± 229	25.92 ± 1.36	1.12	B
306B	⋯	83.81600000	−5.38958333	21.15 ± 0.54	⋯	⋯	3473 ± 35	40.15 ± 2.53	0.64	B
313	198	83.82279167	−5.38919444	26.35 ± 0.19	1.07 ± 0.86	3.21 ± 0.20	4206 ± 235	10.04 ± 0.93	0.64	E2
322	⋯	83.81787500	−5.38794444	24.90 ± 0.27	⋯	⋯	3414 ± 18	37.47 ± 0.92	0.07
324	226	83.81183333	−5.38777778	30.77 ± 0.34	0.70 ± 0.55	−1.87 ± 0.04	3449 ± 36	5.12 ± 1.21	0.01
331	121	83.82604167	−5.38769444	31.83 ± 0.69	1.26 ± 0.24	1.03 ± 0.04	4034 ± 452	18.92 ± 1.17	0.64
332A	183	83.82425000	−5.38766667	25.43 ± 0.32	−1.61 ± 0.17	0.60 ± 0.38	4012 ± 37	16.48 ± 0.43	0.01	BC
370	227	83.81616667	−5.38388889	31.24 ± 0.17	0.01 ± 0.74	−0.46 ± 0.09	3802 ± 45	14.43 ± 0.34	0.66
375	560	83.82075000	−5.38361111	33.19 ± 0.45	0.78 ± 0.13	1.15 ± 0.04	4225 ± 244	26.95 ± 0.70	0.41
386	⋯	83.81362500	−5.38241667	27.44 ± 0.18	⋯	⋯	3673 ± 13	11.20 ± 0.10	0.11
388	532	83.82316667	−5.38244444	30.59 ± 0.46	0.65 ± 0.23	0.35 ± 0.04	4279 ± 142	17.83 ± 1.46	0.52
389	⋯	83.81516667	−5.38233333	30.86 ± 0.15	⋯	⋯	3566 ± 34	26.89 ± 0.46	0.68
408	143	83.83008333	−5.38075000	23.52 ± 0.94	−1.51 ± 0.71	−1.07 ± 0.10	3660 ± 105	50.79 ± 1.85	0.02
410	62	83.82133333	−5.38058333	21.42 ± 1.91	−0.85 ± 0.29	0.32 ± 0.12	3910 ± 94	57.77 ± 2.25	0.03
412	151	83.81808333	−5.38030556	32.93 ± 0.32	1.21 ± 0.20	0.51 ± 0.27	3563 ± 28	31.96 ± 1.61	0.80
413	35	83.81454167	−5.38016667	25.51 ± 0.41	−0.43 ± 0.04	−0.27 ± 0.74	3643 ± 44	28.79 ± 0.55	0.05
420	154	83.81600000	−5.37941667	25.91 ± 0.24	−0.38 ± 0.31	−1.82 ± 0.25	3510 ± 35	24.16 ± 0.67	0.00
425	⋯	83.82479167	−5.37930556	24.94 ± 0.25	⋯	⋯	3576 ± 22	21.83 ± 0.32	0.00
431	700	83.81216667	−5.37752778	33.49 ± 0.31	1.68 ± 0.60	0.05 ± 0.94	3603 ± 22	18.80 ± 0.97	0.65
436	⋯	83.82658333	−5.37708333	32.93 ± 0.08	⋯	⋯	3658 ± 347	16.30 ± 0.90	0.11
440	⋯	83.82233333	−5.37661111	28.52 ± 0.48	⋯	⋯	3715 ± 88	23.57 ± 1.07	1.37
450	⋯	83.82087500	−5.37586111	26.23 ± 0.24	⋯	⋯	3674 ± 54	12.14 ± 0.74	0.02
504	64	83.81395833	−5.37100000	27.48 ± 0.78	−0.29 ± 0.15	−0.32 ± 0.02	3969 ± 295	30.65 ± 1.86	3.79
522A	110	83.81787500	−5.36955556	24.75 ± 0.47	0.03 ± 0.10	−1.97 ± 0.14	3441 ± 24	33.48 ± 1.98	0.31	B
522B	642	83.81787500	−5.36955556	24.54 ± 0.28	−0.69 ± 0.11	0.48 ± 0.33	3170 ± 39	21.94 ± 1.20	0.80	B
546	567	83.81250000	−5.36666667	21.51 ± 0.18	2.25 ± 0.09	0.96 ± 0.11	3793 ± 56	15.32 ± 1.18	0.27	V
703	137	83.80750000	−5.36863889	33.39 ± 0.06	−1.51 ± 0.90	−1.03 ± 0.52	3584 ± 27	9.14 ± 0.41	1.24
713	⋯	83.82795833	−5.38247222	22.95 ± 0.31	⋯	⋯	4988 ± 192	6.75 ± 1.54	1.60

Notes.

^a ID number from Kim et al. (2019). ^b Reported uncertainties also include the 0.058 km s⁻¹ systematic uncertainty between calibration frames and the 0.5% variation found from differing $\mathrm{log}g$ and T_eff. ^c Continuum veiling causes extreme degeneracies with T_eff. Caution should be taken when using derived temperatures with high veiling. ^d Order 33 veiling parameter. ^e In this column, B = double star, previously known in the literature (Hillenbrand 1997; Robberto et al. 2013); BC = new binary candidates, previously reported as single in the literature; E1 = escape group 1; E2 = escape group 2; V = RV variable source.

A machine-readable version of the table is available.

Download table as:  DataTypeset image

Using the framework outlined in Section 5.2, we computed intrinsic velocity dispersions (IVDs) after removing all sources in escape groups 1 and 3. The values from our model were

$$\begin{eqnarray*}\begin{array}{rcl}{\bar{v}}_{\alpha } & = & 0.01\pm 0.16\,\mathrm{km}\,{{\rm{s}}}^{-1},\\ {\bar{v}}_{\delta } & = & -0.09\pm 0.18\,\mathrm{km}\,{{\rm{s}}}^{-1},\\ {\bar{v}}_{r} & = & {27.45}_{-0.22}^{+0.21}\,\mathrm{km}\,{{\rm{s}}}^{-1},\\ {\sigma }_{{v}_{\alpha }} & = & 1.64\pm 0.12\,\mathrm{km}\,{{\rm{s}}}^{-1},\\ {\sigma }_{{v}_{\delta }} & = & 2.03\pm 0.13\,\mathrm{km}\,{{\rm{s}}}^{-1},\\ {\sigma }_{{v}_{r}} & = & {2.56}_{-0.17}^{+0.16}\,\mathrm{km}\,{{\rm{s}}}^{-1},\\ {\rho }_{1} & = & 0.06\pm 0.09,\\ {\rho }_{2} & = & -0.11\pm 0.09,\\ {\rho }_{3} & = & -0.05\pm 0.09.\end{array}\end{eqnarray*}$$

The correlation coefficients are consistent with zero, indicating little to no correlation between velocity components.

Our computed velocity dispersions are roughly consistent with other studies in both the tangential (plane of the sky) direction, (${\sigma }_{{v}_{\alpha }}$, ${\sigma }_{{v}_{\delta }}$) = (1.63 ± 0.04, 2.20 ± 0.06) km s⁻¹ (K19), (${\sigma }_{{v}_{\alpha }}$, ${\sigma }_{{v}_{\delta }}$) = (1.85 ± 0.04, 2.45 ± 0.06) km s⁻¹ (Platais et al. 2020), and (${\sigma }_{{v}_{\alpha }}$, ${\sigma }_{{v}_{\delta }}$) = (1.79 ± 0.12, 2.32 ± 0.10) km s⁻¹ (Jones & Walker 1988), and in the line-of-sight direction, ${\sigma }_{{v}_{r}}\approx 3.1$ km s⁻¹ (Fűrész et al. 2008), ${\sigma }_{{v}_{r}}\approx 2.1$–2.4 km s⁻¹ (Tobin et al. 2009), and ${\sigma }_{{v}_{r}}\approx 2.2$ km s⁻¹ (Da Rio et al. 2017). However, our derived value of ${\sigma }_{{v}_{r}}={2.56}_{-0.17}^{+0.16}$ km s⁻¹ is slightly higher than the value measured in Tobin et al. (2009) and Da Rio et al. (2017). This is possibly due to our closer proximity to the ONC core than previous studies or potentially the target selection of fainter, redder (lower-mass) sources. This may also indicate kinematic structure along the line-of-sight direction, similar to the velocity elongation seen in the north–south direction along the filament (e.g., Jones & Walker 1988). It is also possible that this component potentially suffers from the impact of unresolved binaries, which we evaluate in the next section.

We identified one new binary candidate within our sample ([HC 2000] 332A), along with six other targets that are components of apparent doubles (candidate binaries). At least one source ([HC 2000] 546) shows significant RV variability between the APOGEE and NIRSPAO measurements, although no secondary component was detected in the AO imaging. There is only one epoch of APOGEE data and one epoch of NIRSPAO data for [HC 2000] 546; therefore, it is not possible to determine the orbital parameters of this potential binary without additional observations. Two of our sources have multiepoch NIRSPAO observations ([HC 2000] 229 and [HC 2000] 217); however, the RVs computed at each epoch are consistent with one another within the errors.

5.4.1. Unresolved Binarity

5.4.2. Simulating the Effect of Binaries

To determine the effects of unresolved binarity on the IVD, we performed a test similar to Da Rio et al. (2017) and Karnath et al. (2019). First, we generate a random intrinsic dispersion drawn from a uniform sample between 1 and 4 km s⁻¹. Next, we convolve this dispersion with the measurement error distribution. We use our observed measurement error distribution from the APOGEE+NIRSPEC data set and create an inverse cumulative distribution function, which we randomly sample to build the error distribution. Lastly, we convolve this distribution with a velocity distribution from a set of synthetic binaries and compare the final distribution to our observed distribution.

To generate our synthetic systems, we first generate a distribution of stars with masses between 0.1 and 1 M_☉. We apply a binary fraction (discussed below) to our sample and use the mass ratio distribution from Reggiani & Meyer (2013) to determine component masses within each system. Next, for binary systems, we apply the eccentricity distribution from Duchêne & Kraus (2013) and the period distribution from Raghavan et al.(2010).

For our simulations, we generate 135 systems (the same number of sources in our 3D kinematics sample). These synthetic systems are created using the velbin package (Cottaar & Hénault-Brunet 2014; Foster et al. 2015). This sampling assumes random orbits with random orientations and provides a 1D RV distribution.

This distribution is then compared to our observed APOGEE+NIRSPEC RV distribution, requiring that the standard deviation of the synthetic distribution be within 2σ of the observed distribution's standard deviation. An example of one randomly drawn distribution compared to our observed distribution is shown in Figure 11. The binary contribution to the resulting distribution is typically far less than the intrinsic dispersion and measurement error distribution, a result also noted by Da Rio et al. (2017).

Zoom In Zoom Out Reset image size

We kept the first 10⁵ models that fit our similarity criteria stated above and generated a distribution of IVDs that passed the similarity criteria of our observed RV distribution. We performed this test for binary fractions of 0%, 25%, 50%, 75%, and 100%. Figure 12 shows our simulation results compared to our IVD determined in Section 5.4. Figure 12 illustrates that our measured IVD is consistent with any level of binarity; however, for the line-of-sight component to be consistent with the tangential components, the ONC would need a binary fraction ≳75%, which is inconsistent with observations and modeling results. Therefore, the larger measured ${\sigma }_{{v}_{r}}$ is likely due to formation/evolution rather than binary effects, and the apparent elongation in the line-of-sight component appears to be real rather than a systematic.

Zoom In Zoom Out Reset image size

A number of studies have examined whether the ONC is virialized (e.g., Jones & Walker 1988; Da Rio et al. 2014, 2017; Kim et al. 2019; Kuhn et al. 2019). Da Rio et al. (2014) estimated a 1D mean velocity dispersion of σ_v,1D ≈ 1.73 km s⁻¹ if the ONC is in virial equilibrium. To test for a virialized state, we adopt a methodology similar to Kim et al. (2019), computing the 1D velocity dispersion as a function of radial distance outward.

To balance the small numbers in our sample, we chose radial bins of size 1'. For each bin, we computed the velocity dispersion using the methods outlined in Section 5.2. The values of σ_α and σ_δ were then used to compute the 1D velocity dispersion, i.e., ${\sigma }_{v,1{\rm{D}}}^{2}=({\sigma }_{{v}_{\alpha }}^{2}+{\sigma }_{{v}_{\delta }}^{2})/2$. Our computed velocity dispersions are listed in Table 4 and plotted in Figure 13. We also show the 1D velocity dispersion predicted for models of virial equilibrium using the estimated combined gas and stellar mass from Da Rio et al. (2014, solid line). We assigned a 30% mass uncertainty, similar to Kim et al. (2019), which is illustrated with dotted lines. It should be noted that the model from Da Rio et al. (2014) assumes a spherical potential; however, there is evidence that the potential is closer to an elongated spheroid (e.g., Hillenbrand & Hartmann 1998; Carpenter 2000; Kuhn et al. 2014; Kuznetsova et al. 2015; Megeath et al. 2016).

Zoom In Zoom Out Reset image size

Table 4. Velocity Dispersions as a Function of Distance

Radii	${\sigma }_{{v}_{\alpha }}$	${\sigma }_{{v}_{\delta }}$	σv,rad	${\sigma }_{v,\tan }$	σv,1D	${\sigma }_{{v}_{r}}$	${\sigma }_{v,1{{\rm{D}}}_{3{\rm{D}}}}$	N
(arcmin)	(km s−1)	(km s−1)	(km s−1)	(km s−1)	(km s−1)	(km s−1)	(km s−1)
0–1	${1.98}_{-0.27}^{+0.27}$	${2.55}_{-0.31}^{+0.31}$	${2.22}_{-0.29}^{+0.29}$	${2.38}_{-0.31}^{+0.31}$	${2.28}_{-0.29}^{+0.21}$	${3.17}_{-0.39}^{+0.39}$	${2.61}_{-0.31}^{+0.25}$	30
1–2	${1.71}_{-0.23}^{+0.23}$	${1.95}_{-0.24}^{+0.24}$	${1.46}_{-0.19}^{+0.19}$	${2.12}_{-0.26}^{+0.26}$	${1.84}_{-0.21}^{+0.16}$	${2.41}_{-0.28}^{+0.28}$	${2.04}_{-0.23}^{+0.18}$	32
2–3	${1.74}_{-0.19}^{+0.19}$	${1.86}_{-0.21}^{+0.21}$	${1.77}_{-0.20}^{+0.20}$	${1.85}_{-0.20}^{+0.20}$	${1.80}_{-0.17}^{+0.14}$	${2.56}_{-0.29}^{+0.29}$	${2.08}_{-0.20}^{+0.17}$	42
3–4	${1.45}_{-0.22}^{+0.22}$	${2.22}_{-0.33}^{+0.33}$	${1.88}_{-0.29}^{+0.29}$	${1.94}_{-0.27}^{+0.27}$	${1.87}_{-0.29}^{+0.21}$	${1.92}_{-0.30}^{+0.30}$	${1.89}_{-0.28}^{+0.21}$	22

Download table as:  ASCIITypeset image

The 1D velocity dispersion computed from PMs is extremely consistent with the results of Kim et al. (2019), which is expected because our subsample originated from their catalog. The 1D velocity dispersions favor a virialized state for the majority of the ONC. The radial distribution of ${\sigma }_{{v}_{r}}$ is similar to that from Da Rio et al. (2017, see Figure 12 therein) for the ONC, although they only showed one bin from R ∼ 0–0.4 pc. However, they also found a dispersion 1σ higher than the virial equilibrium model predicts. We also computed the 1D velocity dispersion using all three velocity components, ${\sigma }_{v,1{{\rm{D}}}_{3{\rm{D}}}}^{2}\,=({\sigma }_{{v}_{\alpha }}^{2}+{\sigma }_{{v}_{\delta }}^{2}+{\sigma }_{{v}_{r}}^{2})/3$, listed in Table 4. This measurement marginalizes over potential differences in a single velocity component. Figure 13 (bottom panel) shows that only the bin closest to the ONC core appears elevated from the virialized model; however, this method may wash out features in a single velocity component that deviate from a virialized state. These results add to the growing evidence that the core of the ONC is not fully virialized. Additionally, we do not find evidence of global expansion similar to the results of Kim et al. (2019). Without accurate distances to individual sources, we are not able to explore whether there is expansion along the line-of-sight velocity component.

The Trapezium sits approximately in the middle of the "integral-shaped filament" (ISF; Bally et al. 1987), a long filament of gas with an approximate "S" shape, and 0.1–0.3 pc in front of the filament (e.g., Baldwin et al. 1991; Wen & O'dell 1995; O'Dell 2018; Abel et al. 2019). There has been an observed elongation in the line-of-sight velocity component, which Stutz & Gould (2016) attributed to interactions with the ISF using APOGEE data. The mechanism put forth by Stutz & Gould (2016) to explain the observed elongation in velocity space is the "slingshot" mechanism, where stars born along the filament could be ejected due to the filament undergoing transverse acceleration while the protostar continually accretes mass. The reason for the ejection is that "when the protostar system becomes sufficiently massive to decouple from the filament, it is released" (Stutz & Gould 2016). Such a mechanism would provide an additional contribution to a larger velocity dispersion. However, the slingshot mechanism is dependent on the direction of the transverse acceleration of the filament (i.e., along the line of sight or tangential on the plane of the sky). As the filament runs north–south in the region of the Trapezium, this could provide the mechanism for the observed anisotropy in the tangential velocity components, as well as the radial component. This could be an important effect, as Stutz & Gould (2016) found that the gravity of the background filament likely dominates the gravity field from the stars.

From the analysis of Stutz & Gould (2016), the region we have analyzed here contains primarily stars versus protostars, determined based on their excess IR emission. However, Stutz & Gould (2016) only considered the radial component of motion and not the tangential components. Additionally, their RV sources were obtained from APOGEE, providing very few sources within the central 0.1 pc of the core region where we observe the highest velocity dispersion along the line of sight. As such, there is additional work to determine how the ISF might work to influence the measured velocity dispersions. Although we do not provide an in-depth theoretical analysis here to compare to observational data, one is warranted.

There is a known kinematic anisotropy in the tangential velocities along the north–south direction (e.g., McNamara 1976; Hillenbrand & Hartmann 1998; Da Rio et al. 2014; Kim et al. 2019). However, the deviation from tangential to radial (i.e., toward the center of the cluster) anisotropy (${\sigma }_{\tan }/{\sigma }_{\mathrm{rad}}-1;$ Bellini et al. 2018) of 0.03 ± 0.04 found by Kim et al. (2019) is consistent with isotropic velocities in the tangential–radial velocity space. It should be noted that throughout this section, we use "radial" (v_rad) to indicate motion on the plane of the sky pointing toward the cluster center, "tangential" (v_tan) to indicate motion on the plane of the sky that is tangential to the previously mentioned radial component, and "line-of-sight" (v_r ) to indicate the RV component of motion.

With our 3D kinematic information, we can now compute the ratio of the tangential dispersion to the line-of-sight dispersion (${\sigma }_{{v}_{1{\rm{D}}}}/{\sigma }_{{v}_{r}}$). Figure 14 shows the velocity ratios for ${\sigma }_{{v}_{\alpha }}/{\sigma }_{{v}_{\delta }}$ (top) and ${\sigma }_{{v}_{1{\rm{D}}}}/{\sigma }_{{v}_{r}}$ (bottom). The ${\sigma }_{{v}_{\alpha }}/{\sigma }_{{v}_{\delta }}$ shows a north–south elongation, which is consistent with previous studies, e.g., ${\sigma }_{{\mu }_{\alpha \cos \delta }}/{\sigma }_{{\mu }_{\delta }}=0.74\pm 0.03$ (Kim et al. 2019) and b/a ≈ 0.7 (Da Rio et al. 2014); see also Jones & Walker (1988) and Kuhn et al. (2019).

Zoom In Zoom Out Reset image size

We also decomposed tangential velocities to radial (v_rad) and tangential (v_tan) coordinates, both on the plane of the sky, through a change of basis where the radial axis points toward the ONC CoM. We note that both of these components, v_rad and v_tan, are strictly on the plane of the sky and do not include any line-of-sight velocity components within their vectors. The ratio of ${\sigma }_{{v}_{\mathrm{rad}}}/{\sigma }_{{v}_{\tan }}$ is shown in Figure 14 (middle), where the majority of points are consistent with isotropic dispersions. We also compare our results to the parent population from Kim et al. (2019; open symbols), which shows that our subsample shows variation from the parent population, possibly due to selection effects. Our 1D tangential–to–line-of-sight velocity dispersions (${\sigma }_{{v}_{1{\rm{D}}}}/{\sigma }_{{v}_{r}}$) show significant deviation from unity. The combined value for our total sample is ${\sigma }_{{v}_{1{\rm{D}}}}/{\sigma }_{{v}_{r}}\,=0.78\pm 0.19$. These measurements may indicate an elongation in the velocities along the line-of-sight direction.

Platais et al. (2020) used a diagram of the angle of the PM vector in polar coordinates with respect to the cluster center to examine how the position vector from the center of the ONC and tangential velocity vector were related for fast-moving sources. The expectation is that runaway stars should have angles between these two vectors close to zero, that is, both the position vector and PM vector point radially away from the center of the ONC. We can use a similar analysis to see if there is preferential motion in the core of the ONC. Figure 15 shows the distribution of vectorial angles for all 135 sources in our sample with three components of motion. Random orbits should have no preferential angle; however, there is structure seen in Figure 15. Specifically, there is a preference for sources to have vectorial angles around 90° and −90°. This corresponds to sources whose motion is tangential to a vector pointing radially away from the center of the ONC (normal vector), possibly indicating that there is a rotational preference for sources in the central ONC. In Figure 15 (bottom), we compare the cumulative distribution of vectorial angles for our 3D sample to the vectorial angles for all PM sources within 4' of the ONC CoM from Kim et al. (2019; 671 sources) and a random uniform distribution. The PM sample follows the expected behavior for a random uniform distribution of vectorial angles, which implies a potential bias within our 3D subsample.

Zoom In Zoom Out Reset image size

To test for similarity between our 3D subsample and the parent PM sample from Kim et al. (2019), we performed the Kolmogorov–Smirnov (K-S) and Anderson–Darling (A-D) tests. The K-S test gave a test statistic = 0.115 with a p-value = 0.085, which indicates that we can reject the null hypothesis that the two samples come from the same distribution only at an 8.5% probability. The A-D test gave a test statistic = 0.802 with a p-value = 0.153, indicating that we cannot reject that these two distributions are significantly different, similar to the result of the K-S test. Neither of these test results are extremely significant, indicating that the resulting preferential motion of our 3D sample may not be statistically significant. However, it is worth noting the potential biases introduced into our 3D subsample versus the Kim et al. (2019) sample.

We also performed a comparison test against a random uniform sample using the same tests mentioned earlier. For our 3D subsample, the K-S test gave a test statistic = 0.126 with a p-value = 0.023, which indicates that we can reject the null hypothesis that the kinematic sample comes from a uniform distribution at a 2.3% probability. The A-D test gave a test statistic = 1.786 with a p-value = 0.064, which indicates that the null hypothesis can be rejected at the 6.4% level. Again, neither of these test results are extremely significant, motivating the need for a larger kinematic sample within the core. Comparing the parent sample from Kim et al. (2019) to the random distribution gives a K-S test statistic = 0.019 with a p-value = 0.974 and an A-D test statistic = −0.892 with a p-value > 0.25 (value capped). Both tests indicate that the parent sample is consistent with a random distribution.

Our sample was selected for sources bright enough to be targeted with NIRSPAO (or APOGEE), likely selecting more ONC sources that are in the foreground portion of the cluster rather than deeply embedded sources, which may indicate motion differences as a function of depth in the ONC. Our target list was also curated from the [HC 2000] catalog, selecting sources preferentially thought to have masses ≲1 M_⊙. It is possible that lower-mass sources are showing different kinematics than the bulk motion within the ONC. Kuhn et al. (2019) attempted to measure the bulk rotation of the ONC using tangential kinematics but found no significant rotational preference.

It is important to note that this analysis is done assuming all stars are at the same distance. The Gaia satellite (Gaia Collaboration et al. 2016) is currently obtaining parallaxes for over 1 billion sources; however, it has been noted in previous studies that Gaia measurements in the ONC are unreliable due to the high level of nebulosity (e.g., Kim et al. 2019). Indeed, the number of Gaia Early Data Release 3 (Gaia Collaboration et al. 2021) sources within 1'' of our sources from this study and consistent within the 2σ combined uncertainty in μ_α and μ_δ is only three sources. Therefore, a full 3D study of ONC kinematics is not yet possible and motivates the need for future facilities to obtain parallax measurements for ONC sources.