The 3D Kinematics of the Orion Nebula Cluster. II. Mass-dependent Kinematics of the Inner Cluster

We present the kinematic analysis of 246 stars within 4′ from the center of Orion Nebula Cluster (ONC), the closest massive star cluster with active star formation across the full mass range, which provides valuable insights in the formation and evolution of star cluster on an individual-star basis. High-precision radial velocities and surface temperatures are retrieved from spectra acquired by the NIRSPEC instrument used with adaptive optics (NIRSPAO) on the Keck II 10 m telescope. A 3D kinematic map is then constructed by combining with the proper motions previously measured by the Hubble Space Telescope Advanced Camera for Surveys/WFPC2/WFC3IR and Keck II NIRC2. The measured root-mean-squared velocity dispersion is 2.26 ± 0.08 km s−1, significantly higher than the virial equilibrium’s requirement of 1.73 km s−1, suggesting that the ONC core is supervirial, consistent with previous findings. Energy equipartition is not detected in the cluster. Most notably, the velocity of each star relative to its neighbors is found to be negatively correlated with stellar mass. Low-mass stars moving faster than their surrounding stars in a supervirial cluster suggests that the initial masses of forming stars may be related to their initial kinematic states. Additionally, a clockwise rotation preference is detected. A weak sign of inverse mass segregation is also identified among stars excluding the Trapezium stars, although it could be a sample bias. Finally, this study reports the discovery of four new candidate spectroscopic binary systems.


INTRODUCTION
Star clusters are the primary sites for a multitude of star formation processes observed throughout the universe (Lada & Lada 2003;Gutermuth et al. 2009).Studying the formation and evolution of star clusters is therefore of crucial importance to constrain star formation theories.The kinematics of star clusters provide valuable insights in the processes of their formation and evolution.Despite contemporary observational efforts, much of the details regarding formation processes in star clusters remain to be unveiled (e.g., Krumholz et al. 2014).A particular challenge has been generating models that successfully explain the formation of lowmass stars (M < 0.3 M ⊙ ).Initial models of the competitive accretion process naturally explained the formation of low-mass stars by invoking a natal cluster with similarly-sized cores in which some protostars were cut off from the reservoir of material through violent dynamical interactions (e.g., Bate et al. 2003).However, competitive accretion models have difficulties in explaining the existence of protoplanetary disks and wide binary systems amongst the lowest mass stars (e.g., Burgasser et al. 2007).More recently, simulations in which lowmass stars form along dense filaments of gas infalling into the forming cluster have more successfully predicted the observed multiplicity and disk properties (e.g., Bonnell et al. 2008;Kainulainen et al. 2017).The low-mass stars formed via fragmentation within filaments have high velocities that prevent them from accreting additional materials from the environment.The signature of this formation process can potentially be observed in very young, non-relaxed clusters as a negative correla-tion between velocity and mass.Therefore, observations that probe the kinematics of stars in young clusters have the potential to shed new light on the origin of low mass stars.
The ONC is an optimal target for the study of the formation and evolution of star clusters, as it is the nearest (389 ± 3 pc; Kounkel et al. 2018) active massive stellar nursery.At about 2 Myr (Hillenbrand 1997;Reggiani et al. 2011), the youth and proximity of the ONC makes it ideal for studying the early formation process of a cluster via kinematic measurements, such as radial velocities and proper motions.
Despite its proximity to us, kinematic observations are plagued by the nebulosity and crowding in the region, especially towards the center of the ONC where the Trapezium, a collection of the brightest stars at the heart of ONC, lies.Fűrész et al. (2008) and Tobin et al. (2009) conducted a large-scale radial velocity (RV) survey of 1215 and 1613 stars in the ONC respectively using the multi-fiber echelle spectrograph at the 6.5-m MMT and Magellan telescopes.Kounkel et al. (2016) presented a reanalysis of the data in Tobin et al. (2009) as well as more recent supplementary observations.The Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2017) spectrograph on the 2.5-m Sloan Digital Sky Survey (SDSS; York et al. 2000) telescope has also acquired near-infrared high-resolution spectroscopic data towards the broader Orion Complex region (Da Rio et al. 2016, 2017;Kounkel et al. 2018).However, observations mentioned above have limited coverage near the Trapezium due to its dense, crowded, and highly embedded nature.Most recently, the Gaia Data Release 3 (DR3; Gaia Collaboration et al. 2016Collaboration et al. , 2023) ) provides more complete coverage in the area with more precise astrometric solutions.However, Gaia DR3 lacks the spectroscopic data necessary to infer stellar parameters, such as effective temperature and surface gravity, and as an optical mission is also plagued by nebulosity.Theissen et al. (2022, hereafter T22) present a kinematic analysis of 56 sources within 4 ′ of the ONC center observed by Keck II NIRSPAO and 172 sources observed by APOGEE.The study concludes that the central region of the ONC is not fully virialized by measuring its intrinsic velocity dispersion (IVD).Moreover, the radial IVD is found to be higher than the tangential component as measured by proper motions from the Hubble Space Telescope (HST) + Keck (Kim et al. 2019, hereafter K19).The work presented here expands the sample size of the sources observed by NIRSPAO and presents further kinematic analysis of the region.
In spite of the observations and studies on the ONC over several decades, many questions remain unanswered regarding the kinematics characteristics and the formation process of this young and dense cluster.For instance, the current virial state of the ONC is unclear.Hillenbrand & Hartmann (1998) suggests that some portion of the ONC is already bounded and the cluster will eventually become bound.Velocity dispersion measurements seem to confirm that the ONC is moderately supervirial (e.g., Da Rio et al. 2014, 2017;Kim et al. 2019;Theissen et al. 2022).On the other hand, the most recent study on the ONC advocates that the ONC is likely bound, as previous measurements of the virial parameter may be inflated due to ejections in unstable N-body interactions Kounkel et al. (2022).N-body simulations can rule out neither possibility (Kroupa 2000;Scally et al. 2005).Another fundamental question is the formation mechanism of the cluster as a whole.Kounkel et al. (2022) uses the data from Gaia DR3 and identifies a stellar age gradient as a function of their distance from us, implying that the star formation front is propagating into the star cluster possibly triggered by the shockwave from a supernova in the past.Evidence suggests that the ONC could form in the oscillating integral-shaped filament (ISF), a filamentary gas structure associated with the ONC (Bally et al. 1987), and recently ejected from the ISF (Stutz & Gould 2016;Stutz 2018;Matus Carrillo et al. 2023).Such a mechanism of producing protostars is referred to as slingshot scenario.Moreover, a net infall of young stars along the ISF towards the center is discovered (Kounkel et al. 2022), consistent with the gravitational fragmentation star formation mechanism (Bonnell et al. 2008).
In this work, we use the W. M. Keck Observatory to acquire near-infrared (NIR) high-resolution spectra of 91 sources at the central ONC region, 24 of which are newly observed after T22.Combined with an updated analysis of 172 stars observed by APOGEE, astrometric measurements from Gaia DR3, and proper motions from the HST + Keck (K19), a more thorough analysis of the ONC is made possible in this work, pushing the boundary of our understanding of star cluster formation and evolution.
In Section 2, we introduce the new data from the Keck Observatory and the data reduction processes.Section 3 presents the results of three-dimensional kinematics, including radial velocity modeled from the spectrum and proper motions measured previously, and the derivation of stellar masses.In Section 4 we analyze the massdependent kinematics of the ONC core, including the virial state, velocity-mass relation, effective temperature offset between NIRSPAO and APOGEE, and preferred proper motion direction.In Section 5, we report the identification of candidate single-lined spectroscopic binaries (SB1) Parenago 1837, V* V1337 Ori, V* 1279 ori, and Brun 590.Moreover, we simulate the effect of binaries on the velocity dispersion in the section.The implications of our results on kinematic structure, star formation in the ONC are discussed in Section 6. Mass segregation and binarity in the cluster is also explored in the same section.Lastly, Section 7 gives a summary of this study and an outlook of our future observation and research plans.The NIRSPAO sources were chosen based on their inclusion in the proper motion catalog presented in K19, as our goal is to measure the three-dimensional motion of the stars in our sample.The brightness of the sample ranges from 7.4 to 12.9 magnitude in K band.The total size of the sample targeted with new observations of 91 was primarily driven by observing time constraints from an initial selection of 100 targets in the central ONC.The sources were selected to span a range of magnitudes, which roughly corresponds to a range of temperatures and masses.Since our goal is to assess the kinematic behavior of the lowest mass sources, attention was paid to targeting a sufficient number of faint sources to place statistical constraints of their motion.A sky map of the sources targeted in this study are illustrated in Figure 1.Sources observed in this study are pinpointed by magenta circles.
We supplement the sample observed with an updated analysis of 172 stars observed by APOGEE within 4 ′ of the center of the ONC (T22).APOGEE targets are marked in cyan boxes and the RV survey conducted by Tobin et al. (2009) and reanalyzed by Kounkel et al. (2016) are represented as amber triangles in Figure 1.The dashed white circle indicates the 4 ′ -radius from the center.Throughout this work, we adopted the same location used by (Da Rio et al. 2014) as the center of mass of the ONC: α J2000 = 05 h 35 m 16.26 s , δ J2000 = −05 • 23 ′ 16.4 ′′ .All 240 sources in our sample are located within the circle.A zoom-in on this region is shown in the right panel of Figure 1.As can be seen, previous studies do not have extensive coverage of the central re-gion due to its crowded nature and high level of nebulosity.AO-fed observations in the near-infrared greatly aid in the measurement of individual spectra in this region.

Observation
To observe our 91 targets, we utilized NIRSPEC in conjunction with the Keck II facility laser guide star (LGS) adaptive optics (AO) system (McLean et al. 1998(McLean et al. , 2000;;van Dam et al. 2006;Wizinowich et al. 2006;Martin et al. 2018).NIRSPEC is a near-infrared echelle spectrograph on Keck II.The observations were conducted between 2015 to 2022.The number of sources observed with NIRSPEC with AO (NIRSPAO) increased from 56 to 91 compared to the previous study (T22), a ∼ 63% increase.Exposures with NIRSPAO utilize the 0.041 × 2.26" slit in the the NIRSPEC-7 filter and K-new filter before and after the upgrade, covering the wavelength of 1.839-2.630 and 1.907-2.554,respectively.This wavelength regime covers the carbon monoxide (CO) absorption lines around 2.29-2.40µm, which are present in the spectra of low-mass stars.Moreover, the hydrogen Brackett-γ line at 2.166 µm and the Si, Fe, and Ti lines at 2.18-2.19µmare also within the wavelength range which help inferring the stellar parameters of higher-mass stars.The spectrograph splits the incoming starlight into multiple rows so as to fit in the squareshaped detector and each row is referred to as an order.In this work, the wavelength coverage of each order in our setup of the detector offset is 2.044-2.075µm for order 37, 2.100-2.133µm for order 36, 2.162-2.193µm for order 35, 2.224-2.256µm for order 34, 2.291-2.325µm for order 33, and 2.362-2.382µm for order 32.The CO lines fall within the range of orders 32 − 33, while the Si, Ti, and Fe lines are situated in order 35 for sources with higher effective temperatures.Therefore, we primarily use orders 32, 33, and 35 to sample the stellar parameters.The resolution of the spectra is R ∼ 25000 for data collected before 2019 and R ∼ 35, 000 on or after 2019 as a result of the upgrade on Keck.
While some targets were bright enough to serve as a natural guide star (R ≲ 15), the extinction in this region means that most of our sources required LGS.There are sufficient sources in the region to supply the needed R∼18 tip/tilt guide star requirement within 1' of the target (Wizinowich et al. 2006).For the majority of observations, the target was acquired with PA = 0 • .However, in some cases there were two sources close enough together to position on the slit simultaneously.In those instances, we rotated to an appropriate PA to align both stars to the fall on the slit.
For the majority of the targets, we take four exposures by placing the sources in the slit in an upper-lower-lower-upper sequence, or ABBA dither pattern.In a few cases, the number of frames differs from due to either loss of target or interruption of observation.HD37887, a star of spectral type B9.5IV/V, is used as the calibration star at a similar airmass for telluric wavelength adjustment either before or after a science object is observed.A log of all NIRSPAO observations, including the dates and the total time of source, is given in Table 1.

Data Reduction
NIRSPEC Data Reduction Pipeline (NSDRP)1 is a pipeline specifically designed for reducing NIRSPEC spectra and is optimized for point sources.In this work, data reduction is conducted using a modified version of the NSDRP2 .The modification includes spatial rectification using the object trace instead of the order edge traces, spectral rectification and wavelength calibration using etalon lamps, cosmic-ray cleaning of flats, and bad-pixel cleaning (see Hsu et al. 2021b,a, for details) The steps to reduce the data for each source are briefly summarized below.
1. Median combine the flat frames to generate a master flat frame, which is used to find order edges.
2. Run the modified NSDRP pipeline to reduce all the frames after trimming the spectra edges.
3. Perform initial wavelength calibration for each order using etalon or sky lines of the telluric spectra.
The reduced spectra are then forward modeled for stellar parameters, which will be discussed in Section 2.4.

Spectral Forward Modeling
The reduced spectra are coadded and forwardmodeled to derive the stellar parameters.Instead of modeling each individual exposure of the same source as in T22, we coadd the spectra before forward-modeling.Compared to modeling each individual exposure, coaddition helps reduce the white noise, enhancing the signal-to-noise ratio (SNR) of the data, and saving computational resources for spectral forward modeling.The specific steps of coaddition are summarized below.First, the flux of each exposure of the same source are scaled to match the median flux of the frame with the highest signal-to-noise ratio (SNR).It is worth mentioning that scaling does not affect the modeling results.The noise are scaled by the same factor.Next, the fluxes of all the exposures on the same target are averaged, weighted by the inverse square of the corresponding noise on a pixel-wise basis.The noise associated with the coadded spectrum is calculated from the uncertainty propagation equation for weighted averaging.The calculation of the weighted-averaged flux and the noise are illustrated in Equation 1 and Equation 2, respectively.
where f coadd and σ coadd are the coadded flux and noise.f i and σ i are the flux and noise of the i-th frame for a source.
The coadded spectra are then forward modeled using the Spectral Modeling Analysis and RV Tool (SMART3 , Hsu et al. 2021a).We refer the readers to Hsu et al. (2021a,b) and T22 for a detailed description of the modeling procedure.The steps of modeling spectra are briefly outlined below.
The first step is obtaining a precise absolute wavelength solution.A quadratic polynomial provided by the NSDRP is adopted as the initial wavelength solution4 .A more precise wavelength solution with a systematic uncertainty of 0.058 km s −1 is derived by cross-correlating the spectrum of our A-star calibrator, HD 37887, and a high-resolution reference telluric spectrum (Moehler et al. 2014) in an iterative approach (T22).The coefficients of the polynomial are updated in each iteration by fitting the best wavelength shifts for all cross-correlation windows of 100 pixels.The coadded spectrum is calibrated by the telluric frame in order 32, 33, and 35 with the lowest root-mean-square of the residual for the final wavelength solution.
Next, we use the PHOENIX ACES AGSS COND stellar atmospheric models (Husser et al. 2013) to forwardmodel the coadded stellar spectrum via the Markov chain Monte Carlo (MCMC) method (Butler et al. 1996;Blake et al. 2007Blake et al. , 2008Blake et al. , 2010;;Burgasser et al. 2016) realized by the ensemble sampler emcee (Foreman-Mackey et al. 2013).The flux is modeled by the same function of wavelength as in T22.Note that the surface gravity log g can hardly be constrained from the spectral modeling within the observed wavelength range and is therefore fixed to be 4, which is the expected value for young, low-mass stars at the age of the ONC and is consistent with other studies (e.g., Kounkel et al. 2018).To see how a different log g might affect the stellar mass estimation based on modeled effective temperature, we conducted a test on a subset of sources spanning across the temperature and SNR range by vigorously changing log g to 3.5 and 4.5 respectively.The majority of the resulting stellar masses remain within 1σ from the values under the assumption of log g = 4.Only a few sources deviate by 2σ or more due to a trade-off between surface gravity and temperature with the veiling parameter acting as a tuning knob.Overall, the choice of log g = 4 is rationalized as the stellar masses largely remains consistent under a relatively wide range of log g.In addition, metallicity is set to be 0 based on the average value of the ONC (e.g., D'Orazi et al. 2009).
The free parameters we sample and their corresponding limits and initial distribution are summarized in Table 2.The veiling parameter is defined in the same way as in T22.
Each source is sampled with 100 walkers and 300 steps using the KDEMove, discarding the first 200 steps, as the walkers typically converge within the first 100 steps based on the walker plots.We have also verified the consistency of the results by running the MCMC sampler for 500 steps, thereby ensuring convergence.A fine-tuning sampling with the same number of walkers, steps, and prior distributions follows after removing the pixels where the residual deviates from its median value by more than three standard deviations of itself.The masking removes the remaining bad pixels and cosmic rays from the spectrum that were not re-jected by the NSDRP.The final distribution of the last 100 steps of the 100 walkers are considered as the posterior distribution.We take the median of the posterior distribution as the measured value for each parameter, and half the difference between the 16-th and 84-th percentile, or the 1-σ range for a normal distribution, as the associated uncertainty.Heliocentric RVs are corrected for barycentric motion using the astropy function radial velocity correction.
In addition to emcee, we also attempted to use another Bayesian inference tool, PyMultiNest (Buchner et al. 2014), to sample the posterior distribution of the stellar parameters adopting the limits in Table 2 as the priors.The built-in multimodal nested sampling algorithm is expected to have a better performance in sampling multimodal distributions, which helps with disentangling potential degenerate distributions between the effective temperature and veiling parameter when either one of them is high.However, as we will show in Section 3.1, most of the sources have a low veiling parameter, suppressing the degeneracies.Additionally, we noticed a significant increase in the running time as the number of modeled parameters increases.With 9-10 parameters to sample, emcee turns out to be the better option in terms of computation efficiency, which is what we eventually adopted, as emcee can also sample multimodal distributions unbiasedly.
Reanalyzed results of the APOGEE samples within 4 ′ of the center of the ONC were performed by T22.The identical sampling procedure which is used on the NIRSPEC data is applied to the APOGEE H-band data to model stellar parameters including effective temperatures, rotational velocities, and RVs for consistency in methodology.The results conform well with the SDSS/APOGEE results.
For most of the sources, we model order 32 and 33 simultaneously for the stellar parameters.This is where the CO lines reside for low-mass stars.However, this procedure fails for a few high-temperature sources because the spectra in orders 32 and 33 are mostly flat without any notable features.Therefore, we modeled Brackett-γ, silicon, titanium, and iron lines in order 35 for two high-temperature sources, HC2000 291A and HC2000 337.
From top to bottom, Figure 2 shows an example of the spectrum and model in order 32 and order 33 for HC2000 172, and order 35 for HC2000 337 respectively.The top panel in each figure panel shows the observed spectrum as the gray line, the model as the red line, and the model multiplied by the telluric model as the blue line.CO lines in order 32 and 33 are marked on top of the spectra, with their transparency indicating the corresponding lab intensity according to the HITRAN database5 .We implemented a cut of intensity larger than 10 −25 cm mol −1 for order 32 for visualization purposes.The intensity is then normalized from 0.05 to 0.95 as transparency.Bracket-γ, silicon, titanium, and iron lines are indicated as vertical dashed lines in the case where we modeled order 35 in the bottom figure.The bottom panel in each figure shows the residual as the black line and the noise as the shaded area.Typically the residual is less than 5% of the median flux.According to the modeling results, HC2000 172 is a source with T eff = 3882.5 ± 36.7 K, RV = 39.25 ± 0.29 km s −1 , and v sin i = 23.36 ± 0.44 km s −1 .The stellar parameters for HC2000 337 are T eff = 4652.1 ± 284.4 K, v sin i = 20.24± 10.44 km s −1 , and RV = 31.56± 5.83 km s −1 .The uncertainty of order 35 fitting results are still significant even though it is more than 3 times better than the the 32 and 33 fit.Caution should be taken when using modeled parameters in order 35.
Note that the results might be impacted by the fringing in the spectra, which is the primary contributor to the residuals of spectral modeling (Hsu et al. 2021b;T22).But it is unlikely to change most of our results.

Three-Dimensional Velocities
With the stellar parameters derived from spectral modeling and previous measurements of proper motions and parallaxes, a 3D mass-dependent kinematic map of the ONC core can be constructed.
When there are multiple epochs of observations for the same object in the NIRSPAO sources, we take the average of the stellar parameters weighted by the inverse of the square of the associated uncertainties.In the case where there is a match between NIRSPEC and APOGEE sources within 1 ′′ , the radial velocity derived from the NIRSPAO observation is prioritized over the value from APOGEE since NIRSPAO has a higher resolution(25000 or 35000 versus 22500) and to keep consistency with previous works (e.g., T22).
The NIRSPAO and APOGEE sources are then crossmatched with both the proper motion catalog measured by the HST + Keck (K19), and Gaia DR3 (Gaia Collaboration et al. 2016Collaboration et al. , 2023) ) for parallax within a separation of 1 ′′ to construct a 3D kinematic map of the ONC core.We retrieved the Gaia DR3 data using astroquery (Ginsburg et al. 2019).Due to the low quality of as- b The veiling parameter is defined in the same way as in T22.
trometric measurements with Gaia induced by the nebulosity and extinction in the ONC region, we adopted the same generous quality cut of astrometric gof al < 16 and photometric mean g mag < 16 for Gaia DR3 sources in the region selected for cross matching as in K19.The G magnitude of the cross-matched sources ranges from 7.8-16.0.The proper motion measurements of HST + Keck are prioritized over Gaia DR3 for the same concern of astrometric measurement quality in the latter.
To ensure data quality, we applied several constraints on our sample.First, the NIRSPAO and reanalyzed APOGEE sources with a modeled RV uncertainty of no greater than 5 km s −1 are selected.237 sources out of 246 remain after the RV constraint is applied.Furthermore, considering our estimated distance from the ONC of 389 ± 3 pc (Kounkel et al. 2018) and the limited accuracy of Gaia astrometric solutions in the region, a generous distance constraint on our sources is imposed with a minimum distance of 300 pc, a maximum of 500 pc, and a minimum parallax over error parallax over error of 5. 2 additional source is filtered out in this step, leaving 235 sources in total.For the remaining sources, we assume the same distance of 389 ± 3 pc Kounkel et al. (2018) in the following analysis, as only 100 of the sources have adopted Gaia parallax measurements after the aforementioned quality cut and constraints.Moreover the uncertainties in Gaia parallaxes are too large to be accounted for compared to the size of the ONC.The median uncertainty of the parallax after translated into distance is 8.7 pc, whereas the radius of the ONC is only about 3.7 pc.
We are left with a total number of 235 sources after applying both the RV and distance constraints: 85 NIRSPAO sources within 1.52 ′ from the center; 167 APOGEE sources within 4 ′ ; 17 observed with both instruments.
Note that the proper motion measurements of HST + Keck are in the rest frame of the ONC, whereas the Gaia DR3 is in the absolute frame.We transform the Gaia DR3 proper motion into the same reference frame as in HST + Keck by offsetting the former values by the average of their differences in each direction.The offsets are ∆µ α * , ∆µ δ = (1.60,0.08) mas yr −1 .
Figure 3 visualizes the 3D velocities of the sources6 .The proper motions are characterized by the direction and length of the arrows, while RVs are represented by the color code.Sources moving faster away from us with larger RV values are shown in red, and smaller RVs are shown in blue.Discerning the presence of any kinematic structure solely by visual inspection proves challenging.Therefore, a comprehensive analysis will be conducted in Section 4 to explore this further.
Figure 4 illustrates the comparison of the kinematics measurements between values adopted in this study and previous observations.The left panel shows the proper motion comparison between HST + Keck and Gaia DR3.Proper motion in the right ascension and declination are denoted as µ α * and µ δ , and the measurements from Gaia DR3 and HST + Keck are denoted as DR3 and HK in the subscript respectively.The source on the bottom left, 2M05351094-0524486 or Gaia DR3 3017363547934810112, is the furthest away from consistency because it has a high ruwe of 1.6, an indication  of poor astrometric solution from Gaia.In the scope of this work, we only utilize the Gaia values if HST + Keck RV is not available, as the systematics of the later is better understood and well accounted for.A simi-lar discussion of comparison between HST + Keck and Gaia DR2 measurements can be found in K19.The right panel shows the RV comparison of the matched sources between NIRSPAO and APOGEE.The binary candidates are marked in different colors.
The comparison between the forward-modeled parameters of NIRSPAO sources in this work and T22 is illustrated in Figure 5.Most of the sources have consistent effective temperatures and RVs.The median absolute difference in effective temperature is 29 K, and the maximum difference is 641 K, mostly within the errorbars of T22.Parenago 1837 accounts for the largest well-constrained RV difference, which is shown in blue.Compared with T22, most of the new sources in this work have low veiling parameters indicating less dust absorption, which helps disentangle the aforementioned degeneracy between the effective temperatures and veiling parameters, yielding more confidently constrained modeled stellar parameters.

Stellar Mass Derivation
With the effective temperature obtained from the forward modeling, stellar mass can be interpolated from evolutionary models.We assume an identical age of 2±1 Myr for all sources in the scope of this work.The reason is twofold.First, the estimated age of the ONC is well-established at around 2 Myr (Hillenbrand 1997;Reggiani et al. 2011).Second, the stellar age is only used for the interpolation of stellar mass.Even if we allow a relatively large uncertainty of 50% in the stellar age, it does not affect the mass interpolation for low-mass stars, which are the majority of our sample.Figure 6 illustrates the stellar mass as a function of the effective temperature for stars of different ages under the MIST stellar evolutionary model.Degeneracies do not become significant until the effective temperature exceeds 4500 K. Indeed, 94% of our sources have a lower effective temperature than 4500 K. Therefore, it is justified to assume a stellar age of 2 ± 1 Myr for all sources for the purpose of mass interpolation.
Four different stellar evolutionary models are used for the mass interpolation: the MESA Isochrones & Stellar Tracks (MIST; Dotter 2016; Choi et al. 2016), the BHAC15 model (Baraffe et al. 2015), non-magnetic isochrones in Feiden (2016), and the Palla & Stahler (1999) model.The stellar mass is interpolated using the effective temperature under the assumption of a uniform 2 ± 1 Myr age and 0 metallicity.We calculate the uncertainty in the interpolated mass by taking half of the difference between the highest and lowest mass in the model grid when varying the effective temperature and stellar age within their associated uncertainties, respectively.
Figure 7 illustrates the mass comparison using all of the four models as well as the values from Hillenbrand (1997), which uses the evolutionary model by D 'Antona & Mazzitelli (1994).Most of the temperature-based interpolated masses agree with one another.However, there is a discrepancy between masses in Hillenbrand (1997) and our interpolated values.Effective temperature and luminosity are used to determine the stellar mass in their study.The high nebulosity and reddening in the ONC area make reliable determination of the luminosity challenging.The use of more recent stellar evolutionary models and age instead of luminosity as the proxy to infer the stellar masses in this work accounts for the difference.
We present the derived parameters from modeling for all sources in our sample in Table 3, including RV, v sin i, temperature, and mass.

MASS-DEPENDENT KINEMATICS ANALYSIS
OF THE CENTRAL ONC

Virial State and Energy Equipartition
The ONC was previously found to be supervirial (Scally et al. 2005;Da Rio et al. 2014;Kounkel et al. 2018;T22).In this study, we revisit this assertion for the central ONC with a larger and denser sample than before in this region.
Da Rio et al. ( 2014) derived the theoretical velocity dispersion for the ONC based on the stellar and gas density profiles within 3 pc of the Trapezium center.The stellar density is where r is the radius from the center of the ONC, and the gas density is According to the virial theorem, the velocity dispersion σ at a certain radius r is related to the enclosed mass within if in virial equilibrium as in Equation 5GM r 2 = 2σ 2 . (5) Substituting the mass with the integral of the density over the 4 arcmin (about 0.45 pc) radius, the dependence of velocity dispersion on the radius can be explicitly derived as in Equation 6: (6) We measure the velocity dispersion of our sources assuming there is an intrinsic velocity dispersion along with a measurement uncertainty.That is, the velocity of the i-th source in each direction can be parameterized , where σ (α,δ,r) i denotes the intrinsic velocity dispersion and ϵ (α,δ,r) i denotes the measurement uncertainty.Note that in this analysis we excluded sources whose radial velocity deviates more than 3σ from the mean value to avoid the effects of extreme values, which may be caused by the signalto-noise ratio of the data or unresolved binaries.The accepted radial velocity range is 27.57± 13.21 km s −1 , leaving out 5 sources with 2 being Trapezium stars, which are multiple systems.MCMC forward-modeling is adopted to sample the intrinsic velocity dispersion.The same algorithm is used as in T22 for the modeling.We directly report our updated values below.The velocity dispersion in each direction and the onedimensional (1D) velocity dispersion, defined as librium.Our result of 2.26 km/s is clearly larger than this value (by over 6-σ) and indicates that the ONC is supervirial with a virial ratio (kinetic over potential energy) of q = (σ 1D3D /σ equilibrium ) 2 /2 ∼ 0.85, consistent with the value in Da Rio et al. (2014).Therefore, we reconfirm that the ONC center is supervirial.
Figure 8 shows the velocity dispersion as a function of separation from the center of the cluster.To mitigate against any potential bias in the result due to different ways of binning, the sources are binned equally spaced, i.e. with identical bin width, and equally grouped, i.e. with an almost identical number of sources in each bin.In Figure 8, the left column shows the case of equally spaced, while the right column shows equally grouped binning.From top to bottom, each row shows the 1D root-mean-squared velocity dispersion, proper motion component, and radial velocity component, respectively.It can be seen from the middle row in Figure 8 that the proper motion component is consistent with the virial equilibrium model.However, the radial velocity dispersion (bottom row) is higher than the requirement of virial equilibrium.Consequently, the 1D velocity dispersion of the ONC core sits right above the theoretical prediction of virial equilibrium regardless of the binning method, confirming the result that the ONC is not fully virialized (Da Rio et al. 2014; K19; T22).
HC2000 322 83.816458 −5.397194 3150.4 ± 143.9 11.97 ± 2.70  Separation from Trapezium (arcmin) Another dynamical state we can infer from the velocity dispersion is whether energy equipartition has occurred in the cluster.The velocity dispersion should be inversely proportional to the square root of the stellar mass in a cluster where energy equipartition has already occurred via gravitational interactions.Previously, Hillenbrand & Hartmann (1998) did not see evidence of equipartition in the ONC.Here we re-evaluate this conclusion with our newest data.Figure 9 shows the equally grouped 1D velocity dispersion of all three directions σ 1D3D , 1D velocity dispersion in the proper motion directions σ 1Dpm = (σ 2 RA + σ 2 DEC ) /2, and radial direction as a function of the stellar mass.The −1/2 slope is clearly not present.This is conceivable as the relaxation time of the cluster is estimated to be 6.5 Myr (Hillenbrand & Hartmann 1998), much larger than the age of the cluster of about 2 Myr.Therefore, energy equipartition has not yet taken place in the central ONC, in agreement with Hillenbrand & Hartmann (1998).

Velocity-Mass Relation
With kinematic information and mass estimates for sources in the central ONC, we can look for correlations between masses and velocities, which is indicative of whether stars form via filament fragmentation.Figure 10 shows the velocity of each source relative its neighbors (including itself) within a 0.1 pc (or 53 ′′ ) radius versus their masses derived from the four different evolutionary models described in Section 3.2.The radius within which sources are considered as neighbors can be varied and will be further discussed in Section 4.3.The median number of neighbors for each source is 11.The data are shown in blue.We employ the Gaussian kernel density estimation (KDE) to visualize the distribution of the data, and the value of the estimator is colored in blue in the background.The 84-th percentile of the estimator value is marked as the purple curve.As can be seen, the envelope displays a negative trend on its upper edge.To look at this trend more globally, we show the equally-grouped running average in red.Since the trend of velocity versus mass appears roughly linear, we fit a linear relationship v rel = k * m + b to the data using the scipy.stats.linregressfunction (Virtanen et al. 2020), where v rel is the relative velocity, m is the stellar mass, k and b are the slope and intercept, respectively.To determine the best-fit slope and intercept and their associated uncertainties, we resample the relative velocity and the mass of each source from a normal distribution centered at the observed values with standard deviations being the uncertainty of measurement.Only positive values are kept for the linear fitting.We resampled 10 5 times and each time conducted a linear regression to the data.The median of the recorded slope and intercept are chosen as the best fit values, while half the differences between the 16 and 84-th percentile are considered as their associated uncertainties.The values and uncertainties of the slope k and the intercept b are labeled in the legend in each figure.We utilize the measured values to calculate the p-value for the linear fit with the null hypothesis being the slope is zero.As all of the p-values are less than 0.05, we can safely reject the null hypothesis and conclude that the negative correlation is statistically significant.Apart from the linear fit, we also calculate the Pearson correlation coefficient R as a reflection of the degree of their correlation.The best-fit value and uncertainty of R is determined by resampling in the same way as the linear fit.As can be seen from the figure, all four models display a negative correlation.
We utilized 5-fold cross-validation to verify the statistical significance of the negative correlation between relative velocity and stellar mass and ensure that it is not driven by only a few outlier datapoints.The data is randomly partitioned into 5 equal sized groups, commonly referred to as 'folds'.Linear regression is then performed on 4 groups of the data, leaving out a different group each time.The slope of the resulting 5 linear regressions in the case of the MIST model is shown in Figure 11.The uncertainty-weighted average of the slope across the 5 folds is −0.85, and the standard deviation is 0.22, consistent with the result shown in Figure 10 (a).Therefore, by randomly selecting 5 different sets of 80% of the data but still arriving at the same conclusion, we further validated our finding that the relative velocity is negatively correlated with stellar mass.

Effective Temperature Offset Between NIRSPAO and APOGEE
Despite using the same spectral modeling algorithm, there seems to be a systematic discrepancy in the effective temperatures between the 17 cross-matched NIRSPAO and APOGEE sources in our sample.Figure 12 shows the comparison.Effective temperatures derived from APOGEE spectra are higher than NIRSPAO results, with a weighted average offset of 586 K and a maximum difference of 892 K. Several reasons may cause this difference.First, the spectra of APOGEE and NIR-SPEC are in H and K band, respectively.The spectral features that determine the modeled stellar parameters are therefore different between the two sets of observations.For example, the CO lines are more sensitive to low temperature sources, as the intended objective of the NIRSPAO observation is to obser for low-mass stars.Additionally, considering the highly embedded and crowded nature of the region, extinction and reddening are not identical in H and K bands.A future study will simulate the effect of reddening on temperature estimates.Additionally, to enhance the quality of modeling, especially RV, there is an ongoing effort to model the fringing in the spectrum, which is the primary contributor to the residuals.For the purpose of this study, NIRSPAO estimated effective temperatures are prioritized over APOGEE results.
To evaluate the impact of choosing the NIRSPAO temperatures for mass estimates, we offset the NIRSPAO temperature by the weighted-averaged difference of 586 K to simulate the effect on the velocity-mass relation discussed in Section 4.2. Figure 12 shows the slope of the linear fit before and after the offset as a function of the radius within which sources are considered as neighbors when calculating relative velocities.It can be seen that the negative correlation with mass is weaker but still exists after inflating the NIRSPAO temperature to match the APOGEE values.Either before or after the offset, the negative trend between relative velocity and mass is more evident locally, i.e., a smaller threshold of radius for sources to be considered as neighbors.Discretion is advised for the near-zero slope at larger neighboring radius, as the entire area being analyzed is about 0.45 pc in radius.More sources on the periphery would have incomplete neighbors, which could affect the accuracy of the correlation and the underlying significance.Previous studies identified signs of expansion within the ONC (e.g., Kounkel et al. 2022), and a rotational preference in the proper motions (T22).With a combination of HST and Gaia measurements, we are able to re-evaluate the above findings in greater detail.A polar histogram of the angle between the stellar proper motion vector and the separation vector from the ONC center is shown in Figure 13.Specifically, a positive angle represents clockwise rotation about the ONC center, whereas a negative angle represents counter-clockwise rotation.An angle of zero indicates the source is moving radially outward with respect to the ONC center on the plane of the sky, while −180 • means the source is moving towards the center.More sources are in the 0 • than ±180 • bin agrees with the finding that the ONC is experiencing a slight expansion (Kounkel et al. 2022).The peak around 90 • in Figure 13 illustrates that the center of the cluster is undergoing a clockwise rotation.This is consistent with the finding in Strand (1958).T22 identifies a rotational preference in both +90 • and −90 • , or clockwise and counter-clockwise directions simultaneously.The increase in sample size shows that the ONC core is ac- Relative Velocity (km s 1 ) k1 = 0.49 ± 0.22 k2 = 0.89 ± 0.22 k3 = 0.78 ± 0.24 k4 = 1.17 ± 0.23 k5 = 0.95 ± 0.23 Figure 11.5-fold cross-validation of the negative correlation between relative velocity and stellar mass.The slope and its associated uncertainty for each linear regression k1-k5 is labeled in the legend.This test adopted the stellar masses derived from the MIST model.tually experiencing a clockwise rotation, though a larger sample size would help confirm this trend.bin agrees with the finding that the ONC is experiencing a slight expansion (Kounkel et al. 2022).The peak at 90 • here illustrates that the sources have a preference for clockwise rotation on the plane of the sky around the ONC center.
Among the 23 sources that have multiple epochs of RV measurements, 4 of them exhibited strong variability in their radial velocities.We therefore report the discovery of two candidate binary systems, Paranego 1837, Brun 590, V* V1337 Ori, and V* V1279 Ori.

Parenago 1837
Parenago 1837, or HC2000 546, exhibits variation in its RVs measured in three different epochs, first by APOGEE followed by two observations by NIRSPAO.According to Gaia DR3, it has a G magnitude of 13.54, with BP (blue pass) magnitude of 13.23 and RP (red pass) magnitude of 11.95.The derived stellar mass of the primary is 0.52 ± 0.04 M ⊙ according to the MIST model, assuming the primary light dominates the observed spectra.The RVs after barycentric correction are 29.14 ± 0.38km s −1 , 21.34 ± 0.28 km s −1 , and 25.99 ± 0.19 km s −1 measured on on 2019 February 25 th (UT), 2020 Janurary 21 st (UT), and 2021 October 20 th (UT), respectively.The total time span between the first and last observation ∆t = 638 days.Figure 14   To explore the possible properties of the companion, we sample the possible orbital parameters of the system using a Monte Carlo rejection sampler, The Joker (Price- Whelan et al. 2017).The Joker requires specifying the priors including the period limits, the RV semiamplitude K, and the standard deviations of the velocity trend priors.We limit the periods to between 10 days and twice the observation time span ∆t, or 1276 days.Orbital solutions with arbitrarily long period and large RV variations can be obtained with only 3 epochs of observations.Therefore, an upper limit on the period is a reasonable assumption considering the low mass of the object.The semi-amplitude prior is set to 5 km s −1 , slightly larger than the variation of the RVs.The standard deviations of the velocity priors is set to be a relatively large value of 100 km s −1 to allow the sampler to fully explore the orbital parameter space.Additionally, we replace the default prior distribution of eccentricity with a uniform distribution between 0 and 0.9 to evenly explore the parameter space.We generate 10 5 prior samples and 2126 orbital solutions that match the observed RVs remain after the rejection sampling.The semi-amplitude of the primary's RV variation K is related to stellar masses as where M is the primary mass, m is companion mass, T is the orbital period, a is the semi-major axis, I is the inclination, and e is the eccentricity of the orbit (Murray & Correia 2010).The semi-major axis a can .Allowed companion in the semi-major axiscompanion mass parameter space and 2126 sampled systems from The Joker.Each sampled system is shown as a small gray point.The period range between 10 days and twice the observation time span 2∆t or 1276 days are indicated by the blue solid line at the bottom and the amber dashdotted line at the top, respectively.The left boundary is set by the required variation in orbital velocity which are illustrated in green dashed line for circular orbits, and in red for orbits with eccentricity up to 0.9 (see Equation 11).Assuming a companion smaller than the observed primary mass of 0.52 M⊙ sets the limit shown by the dotted black line to the right.Shaded area are the allowed parameter space within which the companion can reside with the assumptions above.The forbidden periods which are integer fraction of the observation time span ∆t are labeled as purple dotted line.
be expressed as the stellar masses and the orbital period T according to the Kepler's third law where µ = G (M + m) is the standard gravitational parameter.Substituting Equation 9 into Equation 8, K can be expressed as a function of stellar masses, eccentricity, and orbital period Therefore, we can solve for the minimum companion mass when I = π/2 given the K, T , and e for each sampled system from Equation 10with the derived primary mass of 0.50 M ⊙ .Figure 15 shows the distribution of sampled systems and the theoretically allowed region in the semi-major axis-companion mass parameter space.The bottom and top limits are set by the assumed period range between 10 days and 1276 days.The left boundary is set by the constraint that the maximum variation in orbital velocity difference when the orbit is observed from an edge-on direction must exceed the observed amplitude of RV change, i.e., where v max is the maximum orbital velocity difference, ∆v obs is the observed RV difference, and v p , v a denotes the orbital velocity at perihelion and aphelion, respectively (Murray & Dermott 1999).Both the cases when e = 0 and e = 0.9 are plotted in Figure 15 in green and red dashed lines.As can be seen, most sampled systems are within the green line, with a small fraction residing between the green and the red line.Assuming a companion smaller than the observed primary mass of 0.50 M ⊙ sets the limit shown by the dotted black line to the right.The gaps visible in the sampled systems and marked by purple dotted lines are the forbidden periods of integer number fractions of the observation time span ∆t, as the first and last measured RVs are not identical.
876 potential orbital fits with a period greater than ∆t/3 or 213 days but less than 2∆t or 1276 days are plotted in gray in Figure 14.Three different modes of orbits can be clearly seen from the figure, corresponding to periods within ranges of ∆t/4 ∼ ∆t/2, ∆t/2 ∼ ∆t, and ∆t ∼ 2∆t, respectively.
Despite the limited epochs of observation, we are able to infer the approximate mass and separation of the companion thanks to the derived primary mass under reasonable assumptions on the orbital period between 10 days and twice the observation time span of 1276 days.Parenago 1837 is a candidate binary system consisting of a primary of 0.52 ± 0.04 M ⊙ and most likely a companion of ∼ 0.03-0.3M ⊙ , with a separation of less than 2 au.Further observation is needed to more robustly constrain its orbit.

V* V1337 Ori
V* V1337 Ori, or HC2000 214, is another binary candidate.The original APOGEE results have 6 vis-its on the object.The difference in RV is as large as 17.87 km s −1 throughout all the visits, which makes it very likely to be a binary system.However, currently we only have one reanalyzed APOGEE result.Single-epoch reanalysis is required to further constrain its orbit.For consistency, we utilize the reanalyzed APOGEE result of 45.52±1.25 km s −1 , followed by NIRSPAO measurement on 2020 January 20 th (UT) of 36.93 ± 1.05 km s −1 .The two measurements differ by more than 7-sigma.The inferred primary mass is 0.52 ± 0.14 M ⊙ .Due to the limited measurements, it would be challenging to do a similar analysis as Parenago 1837.Therefore, reanalysis on individual previous APOGEE visits or future observations are needed to confirm and constrain this binary system.

V* V1279 Ori
V* V1279 Ori, or HC2000 170, is a source of 0.40 ± 0.03 M ⊙ according to the MIST model.It has a RV measurement of 23.7 ± 1.0 km s −1 in the RV survey by Tobin et al. (2009).Our Keck NIRSPAO measures a RV of 32.70 ± 1.62 km s −1 on 2022 Jan 18 th (UT), more than 5-σ different from each other.

Brun 590
Brun 590, or HC2000 172, is another candidate of 0.60 ± 0.06 M ⊙ according to the MIST model with 2 RV measurements.NIRSPAO measures 23.36 ± 0.29 km s −1 on 2022 Jan 20 th (UT).Different interpretations of APOGEE RV is present in the literature.The reanalyzed APOGEE RV from T22 is 29.26 ± 1.23 km s −1 , while Kounkel et al. (2019) gives 19.320 ± 1.182 km s −1 after removing the systematic effect of temperature and epoch-dependent offsets (Cottaar & Hénault-Brunet 2014).Additional validation and observation is required to confirm whether it is a binary system and constrain the orbit.

Velocity Dispersion
The Keck can spatially resolve binaries with separations larger than 25 mas, or about 10 au at the distance of the ONC (Lacour et al. 2011).Closer binary systems can not be resolved.Currently we do not have a measurement of the close binary fraction in the ONC , as most of the published work has focused on visual binaries.A large fraction of optically unresolved binaries will have a profound effect on the IVD.Therefore, we use the velbin package (Cottaar & Hénault-Brunet 2014;Foster et al. 2015) to simulate the effects of unresolved binaries on our IVD measurements following a similar procedure as described in Da Rio et al. ( 2017 The synthetic simulated RV distribution consists of three components: the systematics, the measurement uncertainty, and the binary offset.We briefly describe the steps to reproduce the three components below.First, a random intrinsic velocity dispersion is drawn from a uniform distribution between 1 and 4.5 km s −1 , which roughly covers a symmetric range on both ends lower and higher than the measured radial velocity dispersion of 2.87 ± 0.15 km s −1 .The systematics are the product of the intrinsic velocity dispersion and a standard normal distribution of length 235.The simulated measurement uncertainty is generated by randomly sampling from the cumulative distribution function (CDF) constructed from the observed measurement error distribution and multiplied by another standard normal distribution.Similarly, a mass distribution is generated by sampling from the CDF constructed by the interpolated mass distribution of the 235 sources.An ensemble of stellar binary systems with uniform distribution of mass ratio and eccentricities is then simulated using this mass distribution with velbin, which in turn gives the radial velocity offset given the binary fraction.
The above process is repeated, and the first 10 5 simulations are kept in which the standard deviation of the synthetic distribution lies within 2σ from the standard deviation of the observed distribution, where σ is the uncertainty in the velocity dispersion in the radial component, or 0.15 km s −1 in our case.Note that radial velocity offsets with an absolute value greater than 7 km s −1 are truncated when fitting for the standard deviations to avoid the impact of extreme values.
Figure 16 shows the simulated and observed velocity distributions in one of the simulations with binary fraction set to 50%. Figure 17 shows the distribution of the 10 5 simulated IVDs that satisfy the aforementioned criterion under different imposed binary fractions.The blue violin plot illustrates the distribution of the simulated radial velocity dispersion.IVDs in the radial, right ascension, and declination directions are shown in red dotted line, green dash-dotted line, and purple dashed lines respectively for comparison.The crossing point between the measured radial velocity dispersion and the interpolation of the median of the simulated radial velocity dispersion at each binary fraction is at 62.54%.In other words, a binary fraction ≳ 62.54% is required to account for our measured higher velocity dispersion in the radial direction if the higher value is solely induced by binaries, which is unreasonably high compared to estimates in the literature.As can be seen from Figure 17, the IVD of the radial component is higher than the other two directions regardless of the simulated binary fraction.Therefore, the binarity alone is not sufficient to explain the larger values in the radial component.

Unresolved Binary Mass
Unresolved close binaries are a potential source of systematic effects that could influence our results.Despite a lack of measurement of spectroscopic binaries in the region, we utilized multiplicity surveys collected in Offner et al. (2023) and conducted a quantitative test on how the unresolved close binaries affect the correlation between relative velocity and stellar mass.Since Keck has a spatial resolution of 10 au at the distance of the ONC, we adopted the close binary fraction (CBF) within 10 au for brown dwarfs and main sequence stars to calculate the mass contribution from unresolved binaries.Historically, the distribution of the mass ratio q is approximated by a power-law f q ∝ q γ .We chose the values of γ for binaries with a separation between 1 au and 10 au.The expectation of the mass of the hidden companion is therefore where the values of CBF and q is determined by the observed primary mass according to Offner et al. (2023).Specifically, we utilized the values from the surveys in Winters et al. (2019); Raghavan et al. (2010), andTokovinin (2014) in orders of increasing stellar mass.Note that γ is assumed to be 0 where the value is unavailable.The closest available values are utilized for sources situated in gaps of stellar mass ranges.For overlapping mass ranges, we adopt the values in the survey where the source is closer to the center of its mass range.Consequently, the adopted values are  12for each source gives the expectation of the mass of the unresolved close binary.As a result, the stellar mass increases by 8-14%.Additionally, the negative correlation in Figure 10 persists after accounting for this excess in mass.The unresolved binaries only produces a marginally flatter linear fit but still with a negative slope.For a quantitative comparison, the slopes in Figure 10 becomes −0.75 ± 0.18, −1.05 ± 0.19, −0.71 ± 0.19, and −0.72 ± 0.14 for MIST, BHAC15, Feiden, and Palla models respectively.This is because the stellar mass on the higher-mass end is shifted further right in Figure 10, while the change in mass of the lower-mass source is limited.Therefore, the excess mass stretches the distri-bution of the data points horizontally, but the negative correlation persists.In this analysis, we have reconfirmed the supervirial nature of the central ONC.This is primarily driven by the measurement of higher velocity dispersion in the radial dimension, which is not due to unresolved binaries.We have also, for the first time, identified a tentative negative trend in relative velocity of stars as a function of mass, with lower mass stars having higher velocities than high mass stars.This has potentially strong implications for star formation.
The primary pathway of stellar formation across the full mass range in stellar clusters remains uncertain.Bonnell et al. (2008) conducted hydrodynamical simulations to investigate low-mass star and brown dwarf formation in clusters.They argue that the filament-shaped infalling gas that is accreted onto a star cluster has high densities, allowing low-mass stars and brown dwarfs to form.However, the high velocity and tidal shear within the gas preclude those low-mass objects from accreting significantly from their surroundings any further.Therefore, one observable feature would be lower mass stars having higher velocities relative to their neighbors and vice versa.The baseline assumptions of this simulation are well-matched to the ONC: a young cluster residing within a filament of gas.Consequently, the ONC can serve as a perfect laboratory for this theory, which could provide keys to the origin of the initial mass function (IMF).
Our measurements of a negative correlation between velocity and mass potentially indicates that the initial mass of forming stars may indeed depend on their initial kinematic states, supporting the gravitational fragmentation mechanism (Bonnell et al. 2008).Note that the negative correlation is more significant in the simulation in Bonnell et al. (2008), as the simulated cluster is still in its nascent age when inspected, only about 0.39 Myr after the first stars formed.The negative trend in the ONC is already partially washed away by dynamical relaxation considering its age of 2 Myr.
Apart from Bonnell et al. (2008), the same trend between velocity and stellar mass is also detected in magnetohydrodynamical simulations of star cluster formation (Mathew & Federrath 2021).The velocities of the sink particles are found to be negatively correlated with their masses shortly after the birth of the stars.The disappearance of the correlation is observed as well over time due to dynamical evolution.We suggest that the ONC is undergoing a similar process of losing the currently observed trend between relative velocity and stellar mass.
Note that the RV measurements of the original APOGEE results were previously found to be negatively-correlated with T eff for low temperature sources below 3400 K (Cottaar & Hénault-Brunet 2014;Kounkel et al. 2019), which would introduce a bias in the velocity mass correlation.However, the reanalyzed APOGEE values adopted in this work are not affected by the bias, as we checked the consistency with the RVs in Kounkel et al. (2019) after removing the systematics.The values conforms well, with a median absolute difference of 0.46 km s −1 .
There are interesting implications of our finding that velocity dispersion in the radial component is larger than that of the proper motion.We speculate that this can be attributed to the influence from the integral-shaped filament (ISF), a gas filament associated with the ONC (Bally et al. 1987).The ISF is believed to experience periodic oscillations in radial and on-sky directions (e.g., Stutz & Gould 2016;Stutz 2018;Matus Carrillo et al. 2023).Protostars are ejected from the gas filaments during the oscillations, shutting off their accretion.Such a process of producing protostars is refereed to as the slingshot mechanism.The ONC may be a star cluster that is radially ejected from the ISF towards us, which results in higher velocities and dispersion in radial direction compared to proper motion components.
The observed expansion of the ONC (Kounkel et al. 2022) and the super-virial velocity dispersion can be explained by the slingshot mechanism as well according to Matus Carrillo et al. (2023).After being ejected from the ISF, the decrease in gas within the ONC reduces the gravitational potential, resulting in the expansion and super-virial velocity dispersion.Alternatively, Kounkel et al. (2022) argues that the expansion is driven by the unstable N-body interactions.Additional observations and tests are required to unveil the reasons for expansion.

3D Spatial Kinematics: Parallax Simulation
Previous analysis was conducted under the assumption that all of the sources are located at the exact same distance of 389±3 pc for reasons discussed in Section 3.1.However, this is not the actual case, especially when we are considering the distances between neighboring sources.Seemingly adjacent sources on the plane of the sky might be distant from one another along our line of sight.
To investigate the impact of unknown distances to individual stars, we chose to simulate the parallax for all of our 235 sources following the same distribution of the limited 100 adopted Gaia parallaxes.An inverse cumulative distribution function (CDF) is constructed from the 100 parallax measurements.We then sample 235 simulated parallaxes from the CDF and assign one value to each source.The simulation is repeated 1000 times.
First, we present the simulation results on the correlation between relative velocity and stellar mass.The distance between sources are now updated to incorporate both the projection on the sky and the radial component.Figure 18 shows the slope k of the linear fit and correlation coefficient R between relative velocity and stellar mass as a function of the separation limit within which we consider sources as neighbors.The blue errorbar and the shaded region represent the mean and standard deviation of the 1000 simulated results of the slope k and correlation coefficient R. As can be seen, the same increasing trend as the separation limit increases persists for both parameters, consistent with Figure 12.Therefore, the conclusion remains based on the 1000 parallax simulations: the negative correlation between relative velocity and stellar mass exists and becomes increasingly obvious when we consider the velocity relative to the more immediate neighbors of each source.
Second, we report the velocity dispersion measurements in the parallax simulation.The same procedure of measuring the velocity dispersion is conducted in each simulation as in Section 4.1.Each of the 1000 simulations produces velocity dispersions and their associated uncertainties.To determine the final uncertainty of the simulation, we first add 1000 normal distributions centering at the simulated values with the standard deviation being the uncertainties of each simulation and normalize it with a factor of 1/1000.Then we fit a single normal distribution to the summation, and use the mean and standard deviation as the simulated values and uncertainties.The summation is justified because we are intrinsically assuming a normal posterior distribution of the parameters when we use a normal distribution likelihood function in the MCMC ensemble sampler.The simulated velocity dispersions are listed below: Comparison of velocity dispersions between the fixed distance scenario and the parallax simulation.The blue profile on the left of each column is a normal distribution of the velocity dispersion when every sources is assumed at a fixed distance of 389 ± 3 pc, with the mean and the standard deviation specified as in Equation 7. The red profile on the right of each column shows the normalized summation of the 1000 normal distributions of the velocity dispersions.
We then fit a single normal distribution over the summed distribution, and the dashed lines show the mean of the normal distributions.The σ1D 3D is shaded with a light gray background to visually distinguish from its components in all three directions.
with simulated parallax is only slightly boosted by 0.03 km s −1 from the case of fixed distance, even smaller than the uncertainty.The change is mostly caused by a minute increase in the RA and DEC components, with the RV component remaining exactly the same.This proves that the projection effect cannot account for the higher σ RV than the other 2 directions.Additionally, the simulated virial ratio is q = 0.88, only slightly different from the previous result of 0.85.In summary, the projection on the plane of the sky does not have a significant effect on the velocity dispersion results.

Mass Segregation and Energy Equipartition
Mass segregation refers to the non-random distribution of stars in stellar systems, where more massive stars tend to concentrate toward the center while lower-mass stars are more dispersed in the outer regions.This occurs due to gravitational interactions and the differential response of stars of different masses to these interactions.Understanding mass segregation and its origin is essential for studying the evolution and formation of stellar systems (e.g., Fregeau et al. 2002;Baumgardt et al. 2008).
With stellar masses derived in Section 3.2, the mass segregation in the central ONC can be quantified.Here we adopt the mass segregation ratio (MSR), Λ MSR , defined in Allison et al. (2009).This parameter uses the minimum spanning tree to measure how clustered the most massive stars in a cluster is.If the most massive stars tend to stay closer to each other compared to a random group of stars of the same amount, the ratio is larger than 1, indicating mass segregation.Conversely, if the most massive stars are further away from each other than a random combination of stars, the ratio is less than 1, or inversely mass segregated.
Figure 20 shows the MSR under the MIST model as a function of N MST .Figure 20 shows that the stars apart from the few Trapezium stars display a feature of inverse mass segregation, i.e., massive stars tend to be on the outskirts.This result may not be a true feature and could instead be a result of our sample selection.As mentioned previously, our goal in the NIRSPAO target selection was to obtain good sampling of low-mass objects in the cluster core.The APOGEE sources are not preferentially sampled at lowmass, tending to higher masses because of higher mass sensitivity limits at H-band.The APOGEE sources are also generally on the periphery of the region we are analyzing.Therefore, we may have missed some sources on the lower-mass end in the outskirts, which results in a bias in this result.More complete spectra in the cluster core and outskirts are needed to assess whether this trend is real or a sample artifact.While the most massive Trapezium stars in the cluster exhibit a close distribution around the cluster center, it is important to recognize that mass segregation does not necessarily indicate the occurrence of energy equiparti-tion, as N-body simulation suggests (Parker et al. 2016).Instead, the kinetic energy decreases at the same rate for low and high-mass stars according to the simulation.This is consistent with our observation in that the stars with higher masses still have a slightly higher velocity dispersion, as can be seen from the left panel in Figure 9.

CONCLUSIONS
In this work, we present the 3D kinematic analysis of 235 sources in the central region of the ONC, including 80 sources observed by the Keck NIRSPAO and 167 reanalyzed APOGEE sources, with 17 common sources between the two surveys.High-precision radial velocities and effective temperatures along with other stellar parameters are retrieved from spectral analysis.With the help of previous proper motion measurements from HST and multiple stellar evolutionary models, we construct a 3D kinematic map of the region with interpolated stellar mass.Listed below is the main takeaways for this work.

A negative correlation between the velocity rel-
ative to the neighbors of each source and the stellar mass is identified using four different stellar evolutionary models, consistent with gravitational filament fragmentation simulation results.This suggests that during the star formation processes within infalling gas filaments, the high velocities of the fast-moving primordial stars preclude them from accreting more materials from their surroundings.
4. Neither the velocity dispersion nor the negative correlation is affected by the projection effect by assuming the same distance of 389 ± 3 pc for all sources.We conducted 1000 simulations in which each source is assigned a simulated distance according to the distribution of the 100 adopted Gaia parallax measurements within our sample.The simulated IVD is consistent with the projected scenario, with σ RA , σ DEC , and σ 1D3D slightly increased by less than 1σ.The negative correlation still holds.
5. There is a systematic discrepancy in the effective temperature between NIRSPAO and APOGEE results.The difference could be attributed to different passband used in the two instruments, K and H band respectively.The negative correlation between relative velocity and stellar mass still exists after accounting for the discrepancy by offsetting the NIRSPAO T eff by the weighted-averaged difference.Furthermore, the negative correlation becomes more significant for more localized relative velocity.
6.A clockwise rotational preference in proper motions is identified in the region.
7. The sources are found to be inversely-mass segregated, or massive stars being more scattered on the outskirts, if the Trapezium stars are excluded.This may stem from a selection bias as we are focusing on low-mass stars closer to the ONC center.
8. We report 4 binary candidate systems by observing the change in radial velocity, Parenago 1837, V* V1337 Ori, V*1279 Ori, and Brun 590.We were able to infer the properties of the companion with a mass most likely of ∼ 0.03-0.3M ⊙ and a semi-major axis most likely of less than 2 au. 9. We simulate the effect of binaries have on the IVD in the radial component and conclude that binarity alone is insufficient to explain the higher value in IVD in the radial direction compared to the proper motion directions.Possible explanations are the interaction with the ISF, such as the slingshot mechansim, which states that the ONC is ejected from the ISF due to the oscillation of the gas filament (Stutz & Gould 2016;Stutz 2018;Matus Carrillo et al. 2023).
10.We calculate the expectation of the excess in mass contributed by unresolved close binaries within 10 au.The expectation of the increase in mass is 8-14%.The negative correlation between relative velocity and stellar mass persists, albeit with a marginally flatter linear fit, after accounting for the contribution from unresolved binary mass.
In the future, we plan to construct kinematic maps in an extended area beyond our current 4 ′ radius.This may place a tighter constraint on the current kinematic state of the larger cluster population, probing deeper into its formation history and implications.In addition, our observations are ideal first epochs for unresolved binary star searches, which will provide a true estimate of the very tight binary fraction in this region.

Figure 1 .
Figure 1.Distribution of sources observed previously and in this study on the background of HST ACS R-band image of the ONC (Robberto et al. 2013).Left: A wide view of the central 30 ′ × 30 ′ of the cluster.The 240 sources observed by NIRSPEC and APOGEE within the circle are selected for analysis in this work.The sources observed by NIRSPAO are marked with magenta circles.Sources observed by APOGEE and by the RV survey in Tobin et al. (2009) are marked in cyan boxes and amber triangles respectively.The dashed white circle indicates the 4 ′ radius threshold.Right: A detailed view within the central 4 ′ radius, highlighting the sources considered in this study.The cluster center is labeled as the blue star.Four candidate spectroscopic binary systems identified in this work, namely Parenago 1837, V* V1337 Ori, V* 1279 Ori, and Brun 590, are marked as green plus, yellow cross, amber diamond, and red hexagon, respectively.

Figure 2 .
Figure 2. Examples of observed and modeled spectra in different orders from NIRSPAO along with atomic and molecular lines.Top: HC2000 172 in order 32; Middle: HC2000 172 in order 33; Bottom: HC2000 337 in order 35.The upper panel in each figure shows the observed spectrum and the model.The vertical lines denote the locations of atomic and molecular spectral lines.The gray line is the normalized observed spectrum flux.The model with and without telluric features are represented as blue and red lines, respectively.The lower panel shows the residual as the black line and the noise as the shaded area.The modeled parameters of each order are labeled beside the corresponding panels.

Figure 3 .
Figure 3. Three-dimensional kinematics of the central 4 ′ of the ONC.The proper motions are denoted by the direction and length of the arrows, and the radial velocities are illustrated by the color.A 1 mas yr −1 key to the quiver plot is shown on the top left, or about 1.844 km s −1 assuming a distance of 389 pc.

Figure 4 .
Figure 4. Kinematics comparison with previous measurements.Left: Proper motion comparison between Gaia DR3 in the absolute frame and HST + Keck in the rest frame.Proper motion in the right ascension and declination are denoted as µα * and µ δ , and the measurements from Gaia DR3 and HST + Keck are denoted as DR3 and HK in the subscript respectively.In the scope of this work, we transformed the Gaia DR3 proper motions into the same reference frame as HST + Keck measurements by offsetting the former ones by the average difference between them in both directions.The offsets in right ascension and declination are ∆µα * , ∆µ δ = (1.60,0.08) mas yr −1 .The red dashed line indicates where the two measurements are equal in RA and DEC directions.Right: RV comparison between sources measured with NIRSPAO and APOGEE.The dashed red line indicates the equal line.3 of the 4 candidate binary systems, Parenago 1837, V* V1337 Ori, and Brun 590, are marked in blue dots, amber square, and green diamond with errorbars respectively.

Figure 5 .Figure 6 .
Figure 5.Comparison between the forward-modeled parameters in this work and T22.Left: Comparison of effective temperature.The median absolute difference is 29 K, with a maximum difference of 641 K.The standard deviation of the differences is 134 K. Middle: Comparison of RV.Note that the most different one shown in blue is the identified binary Parenago 1837.Apart from the binary, the median absolute difference in RV is 0.22 km s −1 , and the standard deviation of the difference itself is 0.49 km s −1 .Right: Comparison of the distribution of veiling parameter of order 33.

Figure 8 .
Figure8.Velocity dispersion as a function of separation from the cluster center.The observed velocity dispersion is shown in red, whereas the requirement of viral equilibrium is illustrated as the black line.The dotted line and the shaded area indicate a 30% total mass error on the virial equilibrium model.The sources are binned equally spaced in the left column, i.e., with identical bin width, and equally grouped in the right column, i.e., with an almost identical number of sources in each bin.The number of sources is labeled on top of each bin.Top: 1D velocity dispersion in all directions; Middle: 1D velocity dispersion of the proper motions; Bottom: Velocity dispersion of the radial velocities.

Figure 9 .
Figure9.Velocity dispersion as a function of stellar mass interpolated from the MIST model.From left to right, the three subfigures show the 1D velocity dispersion of all three directions σ1D 3D , the 1D velocity dispersion of the proper motions σ1D pm , and the radial velocity dispersion σ1D RV respectively.The sources are grouped with equal sizes, and the number of sources in each bin is labeled on top of the corresponding bin.

Figure 10 .
Figure 10.The computed velocity relative to the center of mass of the neighbors of each source within a 0.1-pc (or 53 ′′ ) radius versus stellar mass.Four models are used to interpolate the stellar mass: (a) MIST model; (b) BHAC15 model; (c) Feiden model; (d) Palla model.In each subfigure, the data are represented as the blue points and error bars.The value of the kernel density estimator (KDE) is colored in blue in the background to show the distribution of the data.The purple line marks the 84-th percentile of the KDE value.The equally-grouped running average along with its uncertainties are marked and filled in red.The black line shows the best linear fit to the data, with the slope k and intercept b labeled in the legend.p-value and Pearson's correlation coefficient R is labeled in the bottom right of each figure.All four models display both a negative slope of the linear fit and a negative correlation coefficient.

Figure 12 .
Figure12.Slope of linear fit as a function of separation limit of neighbors.We selected 0.1 pc when analyzing the velocity-mass relation in Section 4.2.Here the separation limit within which sources are counted as neighbors is varied to explore how the linear slope changes.The blue line shows the stellar mass under the MIST model with the derived T eff , while the red line shows the MIST stellar mass after offsetting the NIRSPAO T eff by the average difference of 526.25 K between NIRSPAO and APOGEE.The dashed line shows the zero-slope, or no correlation.The negative correlation between relative velocity and mass is more evident locally within smaller radius when calculating the relative velocity.

Figure 13 .
Figure13.Distribution of the angle between the displacement vector from the center of the ONC and the proper motion vector on the plane of the sky.Positive values stand for clockwise rotation about the center of the ONC, and vice versa.0 • stands for expansion and ±180 • means contraction towards the center.More sources are in the 0 • than ±180 • bin agrees with the finding that the ONC is experiencing a slight expansion(Kounkel et al. 2022).The peak at 90 • here illustrates that the sources have a preference for clockwise rotation on the plane of the sky around the ONC center.
shows each measured RV at the time of observation.The blue and red errorbars indicate the results from APOGEE and NIRSPAO, respectively.

Figure 14 .
Figure 14.The RVs after barycentric correction of the identified binary candidate Parenago 1837 measured by APOGEE and NIRSPAO in three epochs.876 potential orbital fits with a period between ∆t/3 (213 days) and 2∆t (1276 days) are shown in gray.Three different modes of orbits can be clearly seen from the figure, corresponding to periods within ranges of ∆t/4 ∼ ∆t/2, ∆t/2 ∼ ∆t, and ∆t ∼ 2∆t, respectively Figure15.Allowed companion in the semi-major axiscompanion mass parameter space and 2126 sampled systems from The Joker.Each sampled system is shown as a small gray point.The period range between 10 days and twice the observation time span 2∆t or 1276 days are indicated by the blue solid line at the bottom and the amber dashdotted line at the top, respectively.The left boundary is set by the required variation in orbital velocity which are illustrated in green dashed line for circular orbits, and in red for orbits with eccentricity up to 0.9 (see Equation11).Assuming a companion smaller than the observed primary mass of 0.52 M⊙ sets the limit shown by the dotted black line to the right.Shaded area are the allowed parameter space within which the companion can reside with the assumptions above.The forbidden periods which are integer fraction of the observation time span ∆t are labeled as purple dotted line.

Figure 16 .Figure 17 .
Figure16.Radial velocity dispersion in a simulation with binary fraction set to 50%.The red dotted line shows the intrinsic velocity dispersion.The gray histogram and curve shows the observed distribution and fitted normal distribution.The synthetic velocity dispersion with contributions from the intrinsic velocity dispersion, the measurement errors, and the binary offset is shown in blue.The dash-dotted green histogram shows the measurement errors, and the pink dashed histogram shows the contribution from the binaries.

1 .
Kinematic Structure of the ONC and Star Formation Implications

Figure 18 .
Figure 18.Parallax simulation results of the slope of linear fit k and correlation coefficient R between relative velocity and stellar mass under MIST model as a function of separation limits within which sources are considered as neighbors.The blue errorbar and the shaded region show the value and associated uncertainty of the two parameters, while the dashed line shows the zero slope, or no correlation.

Figure 20 .
Figure 20.Mass segregation ratio ΛMSR under MIST model as a function of the number of sources chosen to construct the minimum spanning tree NMST.Left: Including the Trapezium stars.Right: Excluding the Trapezium star.Values greater than 1 indicate mass segregation for the corresponding number of most massive stars, and inverse mass segregation for ΛMSR less than 1.The dividing line is shown in a dashed-red line.combined with the high-mass Trapezium stars, while in Figure 20(b) the Trapezium stars are excluded to unveil the mass segregation for other stars.The dividing line of Λ M SR = 1 is plotted as a red dashed line in both figures.It is evident from Figure 20(a) that the 5 most massive stars are very mass segregated, which is expected as the Trapezium stars are clearly spatially clustered in the ONC center.We utilized the masses of the Trapezium stars from the literature (Weigelt et al. 1999; Close et al. 2012 for θ 1 Orionis A, Vitrichenko et al. 2006 for θ 1 Orionis B, Balega et al. (2014) for θ 1 Orionis C, Allen et al. 2017 for θ 1 Orionis D, and Morales-Calderón et al. 2012 for θ 1 Orionis E).Hillenbrand & Hartmann (1998) argues that the high-mass stars are likely prone to form in the ONC center.Surprisingly, Figure 20(b)shows that the stars apart from the few Trapezium stars display a feature of inverse mass segregation, i.e., massive stars tend to be on the outskirts.This result may not be a true feature and could instead be a result of our sample selection.As mentioned previously, our goal in the NIRSPAO target selection was to obtain good sampling of low-mass objects in the cluster core.The APOGEE sources are not preferentially sampled at lowmass, tending to higher masses because of higher mass sensitivity limits at H-band.The APOGEE sources are also generally on the periphery of the region we are analyzing.Therefore, we may have missed some sources on the lower-mass end in the outskirts, which results in a bias in this result.More complete spectra in the cluster core and outskirts are needed to assess whether this trend is real or a sample artifact.While the most massive Trapezium stars in the cluster exhibit a close distribution around the cluster center, it is important to recognize that mass segregation does not necessarily indicate the occurrence of energy equiparti-

Table 1 .
Log of NIRSPAO Observations

Table 2 .
Spectral Forward-Modeling Free Parameters and Their Bounds

Table 3 .
NIRSPAO Forward-Modeling Results Only a portion of the table is shown here.A complete version of this table is available in the online version of this paper.