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

Lingfeng Wei (魏凌枫) Department of Physics, University of California, San Diego,
La Jolla, CA 92093, USA  Center for Astrophysics and Space Sciences,
University of California, San Diego, La Jolla, CA 92093, USA  Christopher A.
Theissen Department of Astronomy and Astrophysics, University of California,
San Diego, La Jolla, CA 92093, USA  Center for Astrophysics and Space
Sciences, University of California, San Diego, La Jolla, CA 92093, USA  Quinn
M. Konopacky Department of Astronomy and Astrophysics, University of
California, San Diego, La Jolla, CA 92093, USA  Center for Astrophysics and
Space Sciences, University of California, San Diego, La Jolla, CA 92093, USA
Jessica R. Lu Department of Astronomy, University of California, Berkeley,
Berkeley, CA 94720-3411, USA  Chih-Chun Hsu Center for Interdisciplinary
Exploration and Research in Astrophysics (CIERA), Northwestern University,
Evanston, IL 60201, USA  Dongwon Kim Riiid, 10F, Parnas Tower, 521, Teheran-
ro, Gangnam-gu, Seoul 06164, Republic of Korea

###### Abstract

We present the kinematic anaylsis of $246$ stars within $4\arcmin$ 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 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 three-dimensional kinematic
map is then constructed by combining with the proper motions previously
measured by the Hubble Space Telescope (HST) ACS/WFPC2/WFC3IR and Keck II
NIRC2. The measured root-mean-squared velocity dispersion is $2.26\pm
0.08~{}\mathrm{km}\,\mathrm{s}^{-1}$, significantly higher than the virial
equilibrium’s requirement of $1.73~{}\mathrm{km}\,\mathrm{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, though it could be a
sample bias. Finally, this study reports the discovery of four new candidate
spectroscopic binary systems.

Star Formation (1569); Stellar kinematics (1608); Star forming regions (1565);
Star Clusters (1567); Radial velocity (1332); Initial mass function(796)

††facilities: HST(ACS/WFPC2/WFC3IR), SDSS(APOGEE), Keck††software:  astropy
(Astropy Collaboration et al., 2013, 2018, 2022), astroquery (Ginsburg et al.,
2019), smart (Hsu et al., 2021a, b), emcee (Foreman-Mackey et al., 2013),
velbin (Cottaar & Hénault-Brunet, 2014; Foster et al., 2015), TheJoker (Price-
Whelan et al., 2017)

##  1 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 low-mass stars ($M<0.3\,M_{\sun}$).
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 low-mass
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 correlation 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\pm 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. 2016, 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\arcmin$ 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 mass-dependent 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.

##  2 Observation and Data Reduction

###  2.1 Sample Selection

Our sample consists of a total of $246$ sources with 2D positions and 3D
velocities. $91$ of them are observed by NIRSPAO on Keck II within $2\arcmin$
of the central Trapezium, and $172$ sources are re-analyzed APOGEE sources
within $4\arcmin$ of the ONC core, with $17$ sources overlapping.

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\arcmin$ 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\arcmin$-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:
$\alpha_{J2000}=05^{h}35^{m}16.26^{s},\delta_{J2000}=-05^{\circ}23\arcmin
16.4\arcsec$. 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 region
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.

###  2.2 Observation

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\arcmin\times 30\arcmin$ 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\arcmin$ radius threshold. Right: A detailed view
within the central $4\arcmin$ 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.

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, 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 $\sim 63\%$ increase. Exposures with NIRSPAO
utilize the $0.041\times 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~{}\mu m$,
which are present in the spectra of low-mass stars. Moreover, the hydrogen
Brackett-$\gamma$ line at $2.166~{}\mu m$ and the Si, Fe, and Ti lines at
$2.18$–$2.19\mu m$ are 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 square-shaped
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~{}\mu m$ for order $37$, $2.100$–$2.133~{}\mu m$ for order
$36$, $2.162$–$2.193~{}\mu m$ for order $35$, $2.224$–$2.256~{}\mu m$ for
order $34$, $2.291$–$2.325~{}\mu m$ for order $33$, and $2.362$–$2.382~{}\mu
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
$\mathrm{R}\sim 25000$ for data collected before 2019 and $\mathrm{R}\sim
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\lesssim 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$\sim$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
$\mathrm{PA}=0^{\circ}$. 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 $4$ due to either loss of
target or interruption of observation. HD$37887$, 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.

Table 1: Log of NIRSPAO Observations HC2000 ID | Date of | No. of | Exp. Time
---|---|---|---
| Obs. (UT) | FramesaaNumber of ABBA frames. | Per Frame (s)
HC2000 322 | 2015 Dec 23 | 4 | 300
HC2000 296 | 2015 Dec 23 | 4 | 1200
HC2000 259 | 2015 Dec 23 | 4 | 90
HC2000 213 | 2015 Dec 23 | 4 | 60
HC2000 306A | 2015 Dec 24 | 4 | 180
HC2000 306B | 2015 Dec 24 | 4 | 180
HC2000 291A | 2015 Dec 24 | 4 | 600
HC2000 291B | 2015 Dec 24 | 4 | 600
HC2000 252 | 2015 Dec 24 | 4 | 300
HC2000 250 | 2015 Dec 24 | 4 | 1200
HC2000 244 | 2015 Dec 24 | 4 | 600
HC2000 261 | 2015 Dec 24 | 4 | 900
HC2000 248 | 2016 Dec 14 | 4 | 600
HC2000 223 | 2016 Dec 14 | 4 | 300
HC2000 219 | 2016 Dec 14 | 4 | 600
HC2000 324 | 2016 Dec 14 | 3 | 1200
HC2000 295 | 2018 Feb 11 | 4 | 450
HC2000 313 | 2018 Feb 11 | 4 | 180
HC2000 332 | 2018 Feb 11 | 4 | 300
HC2000 331 | 2018 Feb 11 | 4 | 450
HC2000 337 | 2018 Feb 11 | 4 | 60
HC2000 375 | 2018 Feb 11 | 4 | 180
HC2000 388 | 2018 Feb 11 | 4 | 120
HC2000 425 | 2018 Feb 12 | 4 | 60
HC2000 713 | 2018 Feb 12 | 4 | 90
HC2000 408 | 2018 Feb 12 | 4 | 450
HC2000 410 | 2018 Feb 12 | 4 | 600
HC2000 436 | 2018 Feb 12 | 4 | 90
HC2000 442 | 2018 Feb 13 | 4 | 450
HC2000 522A | 2019 Jan 12 | 4 | 450
HC2000 522B | 2019 Jan 12 | 4 | 450
HC2000 145 | 2019 Jan 12 | 4 | 600
HC2000 202 | 2019 Jan 12 | 4 | 120
HC2000 188 | 2019 Jan 12 | 4 | 600
HC2000 302 | 2019 Jan 13 | 4 | 450
HC2000 275 | 2019 Jan 13 | 2 | 450
HC2000 245 | 2019 Jan 13 | 4 | 180
HC2000 258 | 2019 Jan 13 | 4 | 180
HC2000 220 | 2019 Jan 13 | 4 | 600
HC2000 370 | 2019 Jan 16 | 4 | 180
HC2000 389 | 2019 Jan 16 | 4 | 120
HC2000 386 | 2019 Jan 16 | 4 | 120
HC2000 398 | 2019 Jan 16 | 4 | 120
HC2000 413 | 2019 Jan 16 | 4 | 180
HC2000 253 | 2019 Jan 16 | 4 | 120
HC2000 288 | 2019 Jan 17 | 4 | 450
HC2000 420 | 2019 Jan 17 | 4 | 450
HC2000 412 | 2019 Jan 17 | 4 | 450
HC2000 282 | 2019 Jan 17 | 4 | 450
HC2000 217 | 2019 Jan 17 | 4 | 180
HC2000 217 | 2020 Jan 18 | 4 | 120
HC2000 229 | 2020 Jan 18 | 4 | 450
HC2000 228 | 2020 Jan 19 | 4 | 120
HC2000 224 | 2020 Jan 19 | 4 | 90
HC2000 135 | 2020 Jan 19 | 4 | 180
HC2000 440 | 2020 Jan 20 | 4 | 120
HC2000 450 | 2020 Jan 20 | 4 | 300
HC2000 277 | 2020 Jan 20 | 4 | 300
HC2000 204 | 2020 Jan 20 | 4 | 300
HC2000 229 | 2020 Jan 20 | 4 | 450
HC2000 214 | 2020 Jan 20 | 4 | 450
HC2000 215 | 2020 Jan 21 | 4 | 300
HC2000 240 | 2020 Jan 21 | 4 | 300
HC2000 546 | 2020 Jan 21 | 4 | 300
HC2000 504 | 2020 Jan 21 | 4 | 120
HC2000 703 | 2020 Jan 21 | 4 | 300
HC2000 431 | 2020 Jan 21 | 4 | 120
HC2000 229 | 2020 Jan 21 | 3 | 450
HC2000 484 | 2021 Feb 01 | 4 | 300
HC2000 476 | 2021 Feb 01 | 6 | 300
HC2000 546 | 2021 Oct 20 | 4 | 300
HC2000 217 | 2021 Oct 20 | 4 | 120
HC2000 277 | 2021 Oct 20 | 2 | 600
HC2000 435 | 2021 Oct 20 | 3 | 600
HC2000 457 | 2022 Jan 18 | 4 | 750
HC2000 479 | 2022 Jan 18 | 4 | 900
HC2000 490 | 2022 Jan 18 | 4 | 300
HC2000 478 | 2022 Jan 18 | 4 | 420
HC2000 456 | 2022 Jan 18 | 4 | 900
HC2000 170 | 2022 Jan 18 | 4 | 120
HC2000 453 | 2022 Jan 19 | 4 | 600
HC2000 438 | 2022 Jan 19 | 4 | 300
HC2000 530 | 2022 Jan 19 | 4 | 900
HC2000 287 | 2022 Jan 19 | 4 | 300
HC2000 171 | 2022 Jan 19 | 4 | 120
HC2000 238 | 2022 Jan 20 | 4 | 600
HC2000 266 | 2022 Jan 20 | 4 | 600
HC2000 247 | 2022 Jan 20 | 4 | 300
HC2000 172 | 2022 Jan 20 | 4 | 300
HC2000 165 | 2022 Jan 20 | 4 | 450
HC2000 177 | 2022 Jan 20 | 4 | 420
HC2000 163 | 2022 Jan 20 | 3 | 420

###  2.3 Data Reduction

NIRSPEC Data Reduction Pipeline (NSDRP)111https://github.com/Keck-
DataReductionPipelines/NIRSPEC-Data-Reduction-Pipeline 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 NSDRP222https://github.com/ctheissen/NIRSPEC-Data-Reduction-Pipeline. 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. 1.

Median combine the flat frames to generate a master flat frame, which is used
to find order edges.

  2. 2.

Run the modified NSDRP pipeline to reduce all the frames after trimming the
spectra edges.

  3. 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.

###  2.4 Spectral Forward Modeling

The reduced spectra are coadded and forward-modeled 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.

$\displaystyle f_{\mathrm{coadd}}$
$\displaystyle=\frac{\sum_{i}\left(f_{i}/\sigma_{i}^{2}\right)}{\sum_{i}\left(1/\sigma_{i}^{2}\right)}~{},$
(1) $\displaystyle\sigma_{\mathrm{coadd}}^{2}$
$\displaystyle=\sum_{i}\left(\frac{\partial f_{\mathrm{coadd}}}{\partial
f_{i}}\cdot\sigma_{i}\right)^{2}=\frac{1}{\sum_{i}\left(1/\sigma_{i}^{2}\right)}~{},$
(2)

where $f_{\mathrm{coadd}}$ and $\sigma_{\mathrm{coadd}}$ are the coadded flux
and noise. $f_{i}$ and $\sigma_{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 (SMART333https://github.com/chihchunhsu/smart, 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
solution444https://www2.keck.hawaii.edu/koa/nsdrp/documents/NSDRP_Software_Design.pdf.
A more precise wavelength solution with a systematic uncertainty of
$0.058~{}\mathrm{km}\,\mathrm{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 forward-model the coadded stellar spectrum via the Markov chain
Monte Carlo (MCMC) method (Butler et al., 1996; Blake et al., 2007, 2008,
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\sigma$ from the values under the
assumption of $\log g=4$. Only a few sources deviate by $2\sigma$ 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 rejected 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$–$\sigma$ 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\arcmin$ 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.

Table 2: Spectral Forward-Modeling Free Parameters and Their Bounds Stellar Parameter | Notation | Limit | Initial Distributionaa$U\left(a,b\right)$ denotes uniform distribution between $a$ and $b$. | Unit
---|---|---|---|---
Effective Temparature | $T_{\mathrm{eff}}$ | $\left(2300,7000\right)$ | $U\left(2300,7000\right)$ | K
Rotational Velocity | $v\sin{i}$ | $\left(0,100\right)$ | $U\left(0,40\right)$ | $\mathrm{km}\,\mathrm{s}^{-1}$
Radial Velocity (RV) | $V_{r}$ | $\left(-100,100\right)$ | $U\left(-100,100\right)$ | $\mathrm{km}\,\mathrm{s}^{-1}$
Airmass | $AM$ | $\left(1,3\right)$ | $U\left(1,3\right)$ |
Precip. Water Vapor | $PWV$ | $\left(0.5,10\right)$ | $U\left(0.5,10\right)$ | mm
Line Spread Function | $\Delta\nu_{\mathrm{inst}}$ | $\left(1,20\right)$ | $U\left(1,10\right)$ | $\mathrm{km}\,\mathrm{s}^{-1}$
Veiling Param.bbThe veiling parameter is defined in the same way as in T22. | $C_{\mathrm{veil}}$ | $\left(0,10^{20}\right)$ | $U\left(0,10^{5}\right)$ |
Noise Factor | $C_{\mathrm{noise}}$ | $\left(1,50\right)$ | $U\left(1,5\right)$ |
Wavelength Offset | $C_{\lambda}$ | $\left(-1,1\right)$ | $U\left(-0.2,0.2\right)$ | Å

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-$\gamma$, silicon, titanium, and iron lines in
order $35$ for two high-temperature sources, HC$2000~{}291$A and
HC$2000~{}337$.

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.

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 database
555https://www.spectralcalc.com/spectral_browser/db_data.php. We implemented a
cut of intensity larger than $10^{-25}~{}\mathrm{cm}\,\mathrm{mol}^{-1}$ for
order $32$ for visualization purposes. The intensity is then normalized from
$0.05$ to $0.95$ as transparency. Bracket-$\gamma$, 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_{\mathrm{eff}}=3882.5\pm
36.7~{}\mathrm{K}$, $\mathrm{RV}=39.25\pm
0.29~{}\mathrm{km}\,\mathrm{s}^{-1}$, and $v\sin{i}=23.36\pm
0.44~{}\mathrm{km}\,\mathrm{s}^{-1}$. The stellar parameters for HC2000 337
are $T_{\mathrm{eff}}=4652.1\pm 284.4~{}\mathrm{K}$, $v\sin{i}=20.24\pm
10.44~{}\mathrm{km}\,\mathrm{s}^{-1}$, and $\mathrm{RV}=31.56\pm
5.83~{}\mathrm{km}\,\mathrm{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.

##  3 Results

###  3.1 Three-Dimensional Velocities

Figure 3: 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 $\mu_{\alpha^{*}}$ and $\mu_{\delta}$, 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
$\left(\widetilde{\Delta\mu_{\alpha^{*}}},~{}\widetilde{\Delta\mu_{\delta}}\right)=\left(1.60,~{}0.08\right)\mathrm{mas}\,\mathrm{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 4:
Three-dimensional kinematics of the central $4\arcmin$ 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~{}\mathrm{mas}\,\mathrm{yr}^{-1}$ key to the quiver plot is shown on the
top left, or about $1.844~{}\mathrm{km}\,\mathrm{s}^{-1}$ assuming a distance
of $389$ pc.

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\arcsec$, 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 cross-matched with both the proper
motion catalog measured by the HST \+ Keck ( K19, ), and Gaia DR3 (Gaia
Collaboration et al., 2016, 2023) for parallax within a separation of
$1\arcsec$ 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 astrometric 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~{}\mathrm{km}\,\mathrm{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\pm 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\pm 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\arcmin$ from the
center; $167$ APOGEE sources within $4\arcmin$; $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
$\left(\widetilde{\Delta\mu_{\alpha^{*}}},~{}\widetilde{\Delta\mu_{\delta}}\right)=\left(1.60,~{}0.08\right)\mathrm{mas}\,\mathrm{yr}^{-1}$.

Figure 3 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 $\mu_{\alpha^{*}}$ and
$\mu_{\delta}$, 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 similar 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.

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~{}\mathrm{km}\,\mathrm{s}^{-1}$, and the standard deviation of the
difference itself is $0.49~{}\mathrm{km}\,\mathrm{s}^{-1}$. Right: Comparison
of the distribution of veiling parameter of order 33.

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.

Figure 4 visualizes the 3-D velocities of the sources666Note that the proper
motions of right ascension are incorrectly labeled in the inverse direction in
Figure 8 in T22. . 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.

###  3.2 Stellar Mass Derivation

Figure 6: Mass-temperature relation for stars of different ages using the MIST
stellar evolutionary model. The red line and shaded area show the model for
stars of age $2\pm 1$ Myr, which is our assumption for stellar age in this
work. Figure 7: Interpolated stellar mass using four different stellar
evolutionary models: MIST, BHAC15, Feiden, and Palla. Also shown are the
values in Hillenbrand (1997) which uses the evolutionary tracks in D’Antona &
Mazzitelli (1994) for stars less than $3M_{\odot}$.

With the effective temperature obtained from the forward modeling, stellar
mass can be interpolated from evolutionary models. We assume an identical age
of $2\pm 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\pm 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\pm 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.

Table 3: NIRSPAO Forward-Modeling Results

HC2000 ID | K19 ID | APOGEE | Gaia DR3 | $\alpha_{J2000}$ | $\delta_{J2000}$ | $T_{\mathrm{eff}}$ | RVaaHeliocentric corrected RV | $\mu_{\alpha}\cos\delta$ | $\mu_{\delta}$ |  $\cdots$ | $M_{\mathrm{MIST}}$
---|---|---|---|---|---|---|---|---|---|---|---
|  |  |  | (deg) | (deg) | (K) | $\left(\mathrm{km}\,\mathrm{s}^{-1}\right)$ | $\left(\mathrm{mas}\,\mathrm{yr}^{-1}\right)$ | $\left(\mathrm{mas}\,\mathrm{yr}^{-1}\right)$ |  $\cdots$ | $\left(M_{\odot}\right)$
HC2000 322 | $\cdots$ |  $\cdots$ | 3017364063323188864 | $83.817875$ | $-5.387944$ | $3373.5\pm 24.8$ | $24.94\pm 0.57$ | $0.48\pm 0.13$ | $-0.06\pm 0.11$ |  $\cdots$ | $0.286\pm 0.021$
HC2000 259 | $\cdots$ |  $\cdots$ | 3017364059023817600 | $83.821125$ | $-5.392778$ | $3434.2\pm 33.7$ | $28.02\pm 0.64$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | $0.321\pm 0.026$
HC2000 213 | $\cdots$ |  $\cdots$ | 3017363955944598016 | $83.816458$ | $-5.397194$ | $3150.4\pm 143.9$ | $11.97\pm 2.70$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | $0.167\pm 0.067$
HC2000 291B | 211 |  $\cdots$ |  $\cdots$ | $83.816000$ | $-5.390444$ | $3171.1\pm 29.0$ | $29.05\pm 0.38$ | $1.10\pm 0.45$ | $-1.63\pm 0.07$ |  $\cdots$ | $0.176\pm 0.015$
HC2000 252 | 118 |  $\cdots$ | 3017364063325263360 | $83.820750$ | $-5.393639$ | $3162.3\pm 112.5$ | $24.73\pm 1.28$ | $2.97\pm 0.37$ | $2.67\pm 1.38$ |  $\cdots$ | $0.172\pm 0.054$
HC2000 250 | 197 |  $\cdots$ | 3017363960253005056 | $83.816250$ | $-5.393889$ | $2977.7\pm 55.1$ | $30.13\pm 0.42$ | $1.70\pm 0.40$ | $-3.85\pm 1.18$ |  $\cdots$ | $0.100\pm 0.012$
HC2000 244 | 180 |  $\cdots$ | 3017363960246048384 | $83.821083$ | $-5.394389$ | $3392.2\pm 17.7$ | $28.13\pm 0.32$ | $-0.02\pm 0.15$ | $0.83\pm 0.40$ |  $\cdots$ | $0.298\pm 0.016$
HC2000 261 | 206 |  $\cdots$ | 3017363960251927936 | $83.814083$ | $-5.392611$ | $3358.7\pm 17.7$ | $26.29\pm 0.20$ | $-0.72\pm 0.24$ | $0.87\pm 0.12$ |  $\cdots$ | $0.277\pm 0.017$
HC2000 248 | 200 |  $\cdots$ |  $\cdots$ | $83.815667$ | $-5.394000$ | $3379.4\pm 40.0$ | $27.06\pm 1.06$ | $1.81\pm 0.43$ | $-2.83\pm 0.20$ |  $\cdots$ | $0.290\pm 0.030$
HC2000 223 | 164 |  $\cdots$ | 3017363960241964032 | $83.814333$ | $-5.395972$ | $3345.9\pm 20.9$ | $27.35\pm 0.29$ | $1.10\pm 0.08$ | $-0.68\pm 0.23$ |  $\cdots$ | $0.270\pm 0.019$
HC2000 219 | 87 |  $\cdots$ | 3017363960251388288 | $83.813167$ | $-5.396306$ | $3480.2\pm 39.8$ | $23.75\pm 0.60$ | $1.87\pm 0.06$ | $-0.65\pm 0.02$ |  $\cdots$ | $0.348\pm 0.031$
HC2000 324 | 226 |  $\cdots$ | 3017364028960744704 | $83.811833$ | $-5.387778$ | $3435.8\pm 29.3$ | $30.70\pm 0.32$ | $0.70\pm 0.55$ | $-1.87\pm 0.04$ |  $\cdots$ | $0.322\pm 0.023$
HC2000 295 | $\cdots$ |  $\cdots$ | 3017364063325260160 | $83.823208$ | $-5.390250$ | $3506.9\pm 65.0$ | $20.69\pm 0.95$ |  $\cdots$ |  $\cdots$ |  $\cdots$ | $0.364\pm 0.047$
HC2000 313 | 198 |  $\cdots$ | 3017364063328048000 | $83.822792$ | $-5.389194$ | $3564.9\pm 74.1$ | $26.21\pm 0.37$ | $1.07\pm 0.86$ | $3.21\pm 0.20$ |  $\cdots$ | $0.398\pm 0.056$
HC2000 332 | 183 | 2M05353124–0523400 | 3017364059023821696 | $83.824250$ | $-5.387667$ | $3997.3\pm 36.1$ | $25.45\pm 0.27$ | $-1.61\pm 0.17$ | $0.60\pm 0.38$ |  $\cdots$ | $0.685\pm 0.073$
HC2000 331 | 121 |  $\cdots$ | 3017364063325254528 | $83.826042$ | $-5.387694$ | $3597.6\pm 103.1$ | $31.43\pm 0.70$ | $1.26\pm 0.24$ | $1.03\pm 0.04$ |  $\cdots$ | $0.418\pm 0.077$
HC2000 337 | $\cdots$ | 2M05353124–0523400 | 3017364059023825280 | $83.827792$ | $-5.387222$ | $4652.1\pm 284.4$ | $31.91\pm 0.28$ | $-0.22\pm 0.04$ | $0.07\pm 0.04$ |  $\cdots$ | $1.252\pm 0.415$
HC2000 375 | 560 | 2M05353124–0523400 | 3017364063330453504 | $83.820750$ | $-5.383611$ | $4040.9\pm 115.3$ | $32.60\pm 0.56$ | $0.78\pm 0.13$ | $1.15\pm 0.04$ |  $\cdots$ | $0.716\pm 0.131$
HC2000 388 | 532 |  $\cdots$ | 3017364063333454336 | $83.823167$ | $-5.382444$ | $4234.5\pm 141.4$ | $30.44\pm 0.44$ | $0.65\pm 0.23$ | $0.35\pm 0.04$ |  $\cdots$ | $0.857\pm 0.163$
HC2000 425 | $\cdots$ |  $\cdots$ | 3017364127743299328 | $83.824792$ | $-5.379306$ | $3546.2\pm 20.3$ | $24.69\pm 0.33$ | $2.37\pm 0.03$ | $2.67\pm 0.03$ |  $\cdots$ | $0.387\pm 0.022$
$\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$ | $\vdots$

Note. — Only a portion of the table is shown here. A complete version of this
table is available in the online version of this paper.

##  4 Mass-Dependent Kinematics Analysis of the Central ONC

###  4.1 Virial State and Energy Equipartition

Figure 8: 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: 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
$\sigma_{\mathrm{1D_{3D}}}$, the 1D velocity dispersion of the proper motions
$\sigma_{\mathrm{1D_{pm}}}$, and the radial velocity dispersion
$\sigma_{\mathrm{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.

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

$\rho_{\mathrm{star}}\simeq
70M_{\odot}\,\mathrm{pc}^{-3}\left(\frac{r}{pc}\right)^{-2.2},$ (3)

where $r$ is the radius from the center of the ONC, and the gas density is

$\rho_{\mathrm{gas}}\simeq 22M_{\odot}\,\mathrm{pc}^{-3}.$ (4)

According to the virial theorem, the velocity dispersion $\sigma$ at a certain
radius $r$ is related to the enclosed mass within if in virial equilibrium as
in Equation 5

$\frac{GM}{r^{2}}=2\sigma^{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:

$\sigma\left(r\right)=\left[\frac{1}{37}\left(\frac{70}{0.8}\left(\frac{r}{pc}\right)^{-0.2}+\frac{22}{3}\left(\frac{r}{pc}\right)^{2}\right)\right]^{1/2}\mathrm{km}\,\mathrm{s}^{-1}~{}.$
(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 as
$v_{(\alpha,\delta,r)_{i}}\pm\sqrt{\sigma_{{(\alpha,\delta,r)}_{i}}^{2}+\epsilon_{{(\alpha,\delta,r)}_{i}}^{2}}$,
where $\sigma_{{(\alpha,\delta,r)}_{i}}$ denotes the intrinsic velocity
dispersion and $\epsilon_{{(\alpha,\delta,r)}_{i}}$ denotes the measurement
uncertainty. Note that in this analysis we excluded sources whose radial
velocity deviates more than $3\sigma$ from the mean value to avoid the effects
of extreme values, which may be caused by the signal-to-noise ratio of the
data or unresolved binaries. The accepted radial velocity range is $27.57\pm
13.21~{}\mathrm{km}\,\mathrm{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 one-dimensional (1D) velocity
dispersion, defined as
$\sigma_{\mathrm{1D_{3D}}}=\sqrt{\left(\sigma_{\mathrm{RA}}^{2}+\sigma_{\mathrm{DEC}}^{2}+\sigma_{\mathrm{RV}}^{2}\right)/3}$,
are listed in Equation 4.1:

$\displaystyle\sigma_{\mathrm{RA}}=1.73\pm
0.09~{}\mathrm{km}\,\mathrm{s}^{-1},$
$\displaystyle\sigma_{\mathrm{DEC}}=2.03\pm
0.11~{}\mathrm{km}\,\mathrm{s}^{-1},$
$\displaystyle\sigma_{\mathrm{RV}}=2.87\pm
0.16~{}\mathrm{km}\,\mathrm{s}^{-1},$
$\displaystyle\sigma_{\mathrm{1D_{3D}}}=2.26\pm
0.08~{}\mathrm{km}\,\mathrm{s}^{-1}.$ (7)

Da Rio et al. (2014) estimate a RMS velocity dispersion in each direction of
$1.73$ km/s if the ONC is in virial equilibrium. Our result of $2.26$ km/s is
clearly larger than this value (by over $6$–$\sigma$) and indicates that the
ONC is supervirial with a virial ratio (kinetic over potential energy) of
$q=\left(\sigma_{\mathrm{1D_{3D}}}/\sigma_{\mathrm{equilibrium}}\right)^{2}/2\sim
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; ; ).

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 $\sigma_{\mathrm{1D_{3D}}}$, 1D velocity dispersion in the
proper motion directions
$\sigma_{\mathrm{1D_{pm}}}=\sqrt{\left(\sigma_{\mathrm{RA}}^{2}+\sigma_{\mathrm{DEC}}^{2}\right)/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).

###  4.2 Velocity-Mass Relation

Figure 10: The computed velocity relative to the center of mass of the
neighbors of each source within a $0.1$-pc (or $53\arcsec$) 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.

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\arcsec$) 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_{\mathrm{rel}}=k*m+b$ to the data using the
scipy.stats.linregress function (Virtanen et al., 2020), where
$v_{\mathrm{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.

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 $k_{1}$–$k_{5}$ is labeled in the legend. This test
adopted the stellar masses derived from the MIST model.

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.

###  4.3 Effective Temperature Offset Between NIRSPAO and APOGEE

Figure 12: 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_{\mathrm{eff}}$, while
the red line shows the MIST stellar mass after offsetting the NIRSPAO
$T_{\mathrm{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.

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 NIRSPEC 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.

###  4.4 Preferred Proper Motion Direction

Figure 13: 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^{\circ}$ stands for expansion and $\pm 180^{\circ}$ means
contraction towards the center. More sources are in the $0^{\circ}$ than $\pm
180^{\circ}$ bin agrees with the finding that the ONC is experiencing a slight
expansion (Kounkel et al., 2022). The peak at $90^{\circ}$ here illustrates
that the sources have a preference for clockwise rotation on the plane of the
sky around the ONC center.

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 counter-clockwise
rotation about the ONC center, whereas a negative angle represents 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^{\circ}$
means the source is moving towards the center. More sources are in the
$0^{\circ}$ than $\pm 180^{\circ}$ bin agrees with the finding that the ONC is
experiencing a slight expansion (Kounkel et al., 2022). The peak around
$90^{\circ}$ 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^{\circ}$ and
$-90^{\circ}$, or clockwise and counter-clockwise directions simultaneously.
The increase in sample size shows that the ONC core is actually experiencing a
clockwise rotation, though a larger sample size would help confirm this trend.

##  5 Binary Candidates and Binary Simulation

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.

###  5.1 Parenago 1837

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 $\Delta t/3$ ($213$ days) and
$2\Delta 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
$\Delta t/4\sim\Delta t/2$, $\Delta t/2\sim\Delta t$, and $\Delta t\sim
2\Delta t$, respectively Figure 15: Allowed companion in the semi-major axis-
companion 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\Delta t$ or $1276$ days are
indicated by the blue solid line at the bottom and the amber dash-dotted 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_{\odot}$ 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 $\Delta t$ are labeled as purple dotted
line.

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\pm 0.04~{}M_{\odot}$ according to the
MIST model, assuming the primary light dominates the observed spectra. The RVs
after barycentric correction are $29.14\pm 0.38\mathrm{km}\,\mathrm{s^{-1}}$,
$21.34\pm 0.28~{}\mathrm{km}\,\mathrm{s^{-1}}$, and $25.99\pm
0.19~{}\mathrm{km}\,\mathrm{s^{-1}}$ measured on on 2019 February 25th (UT),
2020 Janurary 21st (UT), and 2021 October 20th (UT), respectively. The total
time span between the first and last observation $\Delta t=638$ days. Figure
14 shows each measured RV at the time of observation. The blue and red
errorbars indicate the results from APOGEE and NIRSPAO, respectively.

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 semi-amplitude $K$, and the standard
deviations of the velocity trend priors. We limit the periods to between $10$
days and twice the observation time span $\Delta 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~{}\mathrm{km}\,\mathrm{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~{}\mathrm{km}\,\mathrm{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

$K=\frac{m}{M+m}\frac{2\pi a\sin I}{T\sqrt{1-e^{2}}}~{},$ (8)

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 be expressed as the stellar masses and the orbital period $T$ according to
the Kepler’s third law

$a=\left(\frac{\mu}{4\pi^{2}}T^{2}\right)^{1/3}~{},$ (9)

where $\mu=G\left(M+m\right)$ 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

$K=\frac{m}{M+m}\frac{2\pi\sin
I}{\sqrt{1-e^{2}}}\left(\frac{\mu}{4\pi^{2}T}\right)^{1/3}~{}.$ (10)

Therefore, we can solve for the minimum companion mass when $I=\pi/2$ given
the $K$, $T$, and $e$ for each sampled system from Equation 10 with the
derived primary mass of $0.50~{}M_{\odot}$. 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.,

$\begin{split}\Delta v_{\mathrm{max}}&=v_{p}+v_{a}\\\
&=\frac{m}{M+m}\sqrt{\frac{\mu}{a}}\left(\sqrt{\frac{1+e}{1-e}}+\sqrt{\frac{1-e}{1+e}}\right)\\\
&=2\frac{m}{M+m}\sqrt{\frac{\mu}{a\left(1-e^{2}\right)}}\\\ &\geq\Delta
v_{\mathrm{obs}}~{},\end{split}$ (11)

where $v_{\mathrm{max}}$ is the maximum orbital velocity difference, $\Delta
v_{\mathrm{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_{\odot}$ 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
$\Delta t$, as the first and last measured RVs are not identical.

$876$ potential orbital fits with a period greater than $\Delta t/3$ or $213$
days but less than $2\Delta 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 $\Delta t/4\sim\Delta t/2$, $\Delta
t/2\sim\Delta t$, and $\Delta t\sim 2\Delta 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\pm 0.04~{}M_{\odot}$ and most
likely a companion of $\sim 0.03$–$0.3~{}M_{\odot}$, with a separation of less
than $2$ au. Further observation is needed to more robustly constrain its
orbit.

###  5.2 V* V1337 Ori

V* V1337 Ori, or HC2000 214, is another binary candidate. The original APOGEE
results have $6$ visits on the object. The difference in RV is as large as
$17.87~{}\mathrm{km}\,\mathrm{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\pm 1.25~{}\mathrm{km}\,\mathrm{s}^{-1}$, followed by NIRSPAO
measurement on 2020 January 20th (UT) of $36.93\pm
1.05~{}\mathrm{km}\,\mathrm{s}^{-1}$. The two measurements differ by more than
$7$-sigma. The inferred primary mass is $0.52\pm 0.14~{}M_{\odot}$. 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.

###  5.3 V* V1279 Ori

V* V1279 Ori, or HC2000 170, is a source of $0.40\pm 0.03~{}M_{\odot}$
according to the MIST model. It has a RV measurement of $23.7\pm
1.0~{}\mathrm{km}\,\mathrm{s}^{-1}$ in the RV survey by Tobin et al. (2009).
Our Keck NIRSPAO measures a RV of $32.70\pm
1.62~{}\mathrm{km}\,\mathrm{s}^{-1}$ on 2022 Jan 18th (UT), more than
$5$–$\sigma$ different from each other.

###  5.4 Brun 590

Brun 590, or HC2000 172, is another candidate of $0.60\pm 0.06~{}M_{\odot}$
according to the MIST model with $2$ RV measurements. NIRSPAO measures
$23.36\pm 0.29~{}\mathrm{km}\,\mathrm{s}^{-1}$ on 2022 Jan 20th (UT).
Different interpretations of APOGEE RV is present in the literature. The
reanalyzed APOGEE RV from T22 is $29.26\pm
1.23~{}\mathrm{km}\,\mathrm{s}^{-1}$, while Kounkel et al. (2019) gives
$19.320\pm 1.182~{}\mathrm{km}\,\mathrm{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.

###  5.5 Binary Simulation

####  5.5.1 Velocity Dispersion

Figure 16: 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. Figure 17: Binary simulation result of the intrinsic radial velocity
dispersion. The blue violin plot shows the distribution of the $10^{5}$
simulated synthesized radial velocity dispersion under different imposed
binary fractions, with the $0$-th, $16$-th, $50$-th, $84$-th, and $100$-th
percentile marked by the horizontal bars. The measured velocity dispersion in
radial, right ascension, and declination directions are shown in red dotted
line, green dash-dotted line, and purple dashed lines respectively.

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), Karnath et al. (2019), and T22. $235$
systems are generated across the binary fraction range from $0$ to $1$ with a
step size of $0.25$, as is the total number of NIRSPEC and APOGEE sources used
to generate the IVD.

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~{}\mathrm{km}\,\mathrm{s}^{-1}$, which roughly covers a symmetric
range on both ends lower and higher than the measured radial velocity
dispersion of $2.87\pm 0.15~{}\mathrm{km}\,\mathrm{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\sigma$ from the standard deviation of the observed distribution, where
$\sigma$ is the uncertainty in the velocity dispersion in the radial
component, or $0.15~{}\mathrm{km}\,\mathrm{s}^{-1}$ in our case. Note that
radial velocity offsets with an absolute value greater than
$7~{}\mathrm{km}\,\mathrm{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 $\gtrsim~{}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.

####  5.5.2 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}\propto q^{\gamma}$. We chose the values of $\gamma$ for
binaries with a separation between $1$ au and $10$ au. The expectation of the
mass of the hidden companion is therefore

$\mathrm{CBF}\times M\times\frac{\int_{0}^{1}q\cdot
f_{q}dq}{\int_{0}^{1}f_{q}dq}~{},$ (12)

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), and Tokovinin
(2014) in orders of increasing stellar mass. Note that $\gamma$ 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

$\displaystyle\mathrm{CBF}=16\%,$ $\displaystyle\gamma=0$
$\displaystyle\mathrm{if}\quad 0.075<M<0.15~{},$ (13)
$\displaystyle\mathrm{CBF}=14\%,$ $\displaystyle\gamma=0.7$
$\displaystyle\mathrm{if}\quad 0.15<M<0.3~{},$
$\displaystyle\mathrm{CBF}=15\%,$ $\displaystyle\gamma=0.1$
$\displaystyle\mathrm{if}\quad 0.3<M<0.675~{},$
$\displaystyle\mathrm{CBF}=20\%,$ $\displaystyle\gamma=0.2$
$\displaystyle\mathrm{if}\quad 0.675<M<1.0875~{},$
$\displaystyle\mathrm{CBF}=14\%,$ $\displaystyle\gamma=0.4$
$\displaystyle\mathrm{if}\quad 1.0875<M<1.5~{}.$

Substituting the above values back into Equation 12 for 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\pm 0.18$, $-1.05\pm 0.19$, $-0.71\pm 0.19$, and $-0.72\pm 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 distribution of the data points horizontally, but
the negative correlation persists.

##  6 Discussion

###  6.1 Kinematic Structure of the ONC and Star Formation Implications

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_{\mathrm{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~{}\mathrm{km}\,\mathrm{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.

###  6.2 3D Spatial Kinematics: Parallax Simulation

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.

Previous analysis was conducted under the assumption that all of the sources
are located at the exact same distance of $389\pm 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.

Figure 19: 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\pm 3$ pc, with the mean and the
standard deviation specified as in Equation 4.1. 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 $\sigma_{\mathrm{1D_{3D}}}$ is shaded with a
light gray background to visually distinguish from its components in all three
directions.

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:

$\displaystyle\overline{\sigma}_{\mathrm{RA}}=1.78\pm
0.10~{}\mathrm{km}\,\mathrm{s}^{-1},$
$\displaystyle\overline{\sigma}_{\mathrm{DEC}}=2.09\pm
0.11~{}\mathrm{km}\,\mathrm{s}^{-1},$
$\displaystyle\overline{\sigma}_{\mathrm{RV}}=2.87\pm
0.15~{}\mathrm{km}\,\mathrm{s}^{-1},$
$\displaystyle\overline{\sigma}_{\mathrm{1D_{3D}}}=2.29\pm
0.08~{}\mathrm{km}\,\mathrm{s}^{-1}.$ (14)

A comparison of the velocity dispersions between the fixed and simulated
distance scenario is shown in Figure 19. As can be seen, the 1D velocity
dispersion with simulated parallax is only slightly boosted by
$0.03~{}\mathrm{km}\,\mathrm{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
$\sigma_{\mathrm{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.

###  6.3 Mass Segregation and Energy Equipartition

Figure 20: Mass segregation ratio $\Lambda_{\mathrm{MSR}}$ under MIST model as
a function of the number of sources chosen to construct the minimum spanning
tree $N_{\mathrm{MST}}$. 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
$\Lambda_{\mathrm{MSR}}$ less than $1$. The dividing line is shown in a
dashed-red line.

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),
$\Lambda_{\mathrm{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_{\mathrm{MST}}$. Figure 20(a) shows our target sources 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
$\Lambda_{MSR}=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 $\theta^{1}$
Orionis A, Vitrichenko et al. 2006 for $\theta^{1}$ Orionis B, Balega et al.
(2014) for $\theta^{1}$ Orionis C, Allen et al. 2017 for $\theta^{1}$ Orionis
D, and Morales-Calderón et al. 2012 for $\theta^{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 low-mass, 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
equipartition, 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.

##  7 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.

  1. 1.

We model stellar parameters including $T_{\mathrm{eff}}$ and RV using high-
resolution spectra from NIRSPAO ($R\sim 25000$ or $35000$). Stellar masses are
then derived from $T_{\mathrm{eff}}$ and four different stellar evolutionary
models assuming a stellar age of $2\pm 1$ Myr.

  2. 2.

The velocity dispersion in the right ascension, declination, radial direction
are $\sigma_{\mathrm{RA}}=1.73\pm
0.09~{}\mathrm{km}\,\mathrm{s}^{-1},~{}\sigma_{\mathrm{DEC}}=2.03\pm
0.11~{}\mathrm{km}\,\mathrm{s}^{-1},~{}\sigma_{\mathrm{RV}}=2.87\pm
0.16~{}\mathrm{km}\,\mathrm{s}^{-1}$. The 1D velocity dispersion is
$\sigma_{\mathrm{1D_{3D}}}=2.26\pm 0.08~{}\mathrm{km}\,\mathrm{s}^{-1}$,
exceeding the requirement of virial equilibrium of
$1.73~{}\mathrm{km}\,\mathrm{s}^{-1}$ (Da Rio et al., 2014). Therefore, the
ONC is not fully virialized yet, with a virial ratio of $0.85$, consistent
with Da Rio et al. (2014).

  3. 3.

A negative correlation between the velocity relative 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. 4.

Neither the velocity dispersion nor the negative correlation is affected by
the projection effect by assuming the same distance of $389\pm 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 $\sigma_{\mathrm{RA}}$, $\sigma_{\mathrm{DEC}}$,
and $\sigma_{\mathrm{1D_{3D}}}$ slightly increased by less than $1\sigma$. The
negative correlation still holds.

  5. 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_{\mathrm{eff}}$ by
the weighted-averaged difference. Furthermore, the negative correlation
becomes more significant for more localized relative velocity.

  6. 6.

A clockwise rotational preference in proper motions is identified in the
region.

  7. 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. 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 $\sim
0.03$–$0.3~{}M_{\odot}$ and a semi-major axis most likely of less than $2$ au.

  9. 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. 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\arcmin$ 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.

##  8 Acknowledgments

We would like to thank an anonymous referee whose suggestions greatly improved
this manuscript. L.W. thanks the LSSTC Data Science Fellowship Program, which
is funded by LSSTC, NSF Cybertraining Grant #1829740, the Brinson Foundation,
and the Moore Foundation; his participation in the program has benefited this
work. Portions of this work were funded by NSF grant AST-1714816.

The data presented herein were obtained at the W. M. Keck Observatory, which
is operated as a scientific partnership among the California Institute of
Technology, the University of California and the National Aeronautics and
Space Administration. The Observatory was made possible by the generous
financial support of the W. M. Keck Foundation. The authors wish to recognize
and acknowledge the very significant cultural role and reverence that the
summit of Maunakea has always had within the indigenous Hawaiian community. We
are most fortunate to have the opportunity to conduct observations from this
mountain.

The authors thank support astronomers Carlos Alvarez, Greg Doppmann, Jim Lyke,
Luca Rizzi, Hien Tran and telescope operators Joel Aycock, Tony Connors,
Heather Hershley, Julie Renaud-Kim, Gary Puniwai, Arina Rostopchina, Terry
Stickel, and Cynthia Wilburn for their help in obtaining the observations.

Portions of this work were conducted at the University of California, San
Diego, which was built on the unceded territory of the Kumeyaay Nation, whose
people continue to maintain their political sovereignty and cultural
traditions as vital members of the San Diego community.

This work has made use of data from the European Space Agency (ESA) mission
Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing
and Analysis Consortium (DPAC,
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P.
Sloan Foundation, the U.S. Department of Energy Office of Science, and the
Participating Institutions. SDSS acknowledges support and resources from the
Center for High-Performance Computing at the University of Utah. The SDSS web
site is www.sdss4.org.

SDSS is managed by the Astrophysical Research Consortium for the Participating
Institutions of the SDSS Collaboration including the Brazilian Participation
Group, the Carnegie Institution for Science, Carnegie Mellon University,
Center for Astrophysics — Harvard & Smithsonian (CfA), the Chilean
Participation Group, the French Participation Group, Instituto de Astrofísica
de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and
Mathematics of the Universe (IPMU) / University of Tokyo, the Korean
Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut
für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA
Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-
Institut für Extraterrestrische Physik (MPE), National Astronomical
Observatories of China, New Mexico State University, New York University,
University of Notre Dame, Observatório Nacional / MCTI, The Ohio State
University, Pennsylvania State University, Shanghai Astronomical Observatory,
United Kingdom Participation Group, Universidad Nacional Autónoma de México,
University of Arizona, University of Colorado Boulder, University of Oxford,
University of Portsmouth, University of Utah, University of Virginia,
University of Washington, University of Wisconsin, Vanderbilt University, and
Yale University.

## References

  * Allen et al. (2017) Allen, C., Costero, R., Ruelas-Mayorga, A., & Sánchez, L. J. 2017, MNRAS, 466, 4937, doi: 10.1093/mnras/stx076
  * Allison et al. (2009) Allison, R. J., Goodwin, S. P., Parker, R. J., et al. 2009, MNRAS, 395, 1449, doi: 10.1111/j.1365-2966.2009.14508.x
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33, doi: 10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123, doi: 10.3847/1538-3881/aabc4f
  * Astropy Collaboration et al. (2022) Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167, doi: 10.3847/1538-4357/ac7c74
  * Balega et al. (2014) Balega, Y. Y., Chentsov, E. L., Leushin, V. V., Rzaev, A. K., & Weigelt, G. 2014, Astrophysical Bulletin, 69, 46, doi: 10.1134/S1990341314010052
  * Bally et al. (1987) Bally, J., Langer, W. D., Stark, A. A., & Wilson, R. W. 1987, ApJ, 312, L45, doi: 10.1086/184817
  * Baraffe et al. (2015) Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42, doi: 10.1051/0004-6361/201425481
  * Bate et al. (2003) Bate, M. R., Bonnell, I. A., & Bromm, V. 2003, MNRAS, 339, 577, doi: 10.1046/j.1365-8711.2003.06210.x
  * Baumgardt et al. (2008) Baumgardt, H., De Marchi, G., & Kroupa, P. 2008, ApJ, 685, 247, doi: 10.1086/590488
  * Blake et al. (2010) Blake, C. H., Charbonneau, D., & White, R. J. 2010, ApJ, 723, 684, doi: 10.1088/0004-637X/723/1/684
  * Blake et al. (2007) Blake, C. H., Charbonneau, D., White, R. J., Marley, M. S., & Saumon, D. 2007, ApJ, 666, 1198, doi: 10.1086/520124
  * Blake et al. (2008) Blake, C. H., Charbonneau, D., White, R. J., et al. 2008, ApJ, 678, L125, doi: 10.1086/588754
  * Bonnell et al. (2008) Bonnell, I. A., Clark, P., & Bate, M. R. 2008, MNRAS, 389, 1556, doi: 10.1111/j.1365-2966.2008.13679.x
  * Buchner et al. (2014) Buchner, J., Georgakakis, A., Nandra, K., et al. 2014, A&A, 564, A125, doi: 10.1051/0004-6361/201322971
  * Burgasser et al. (2016) Burgasser, A. J., Blake, C. H., Gelino, C. R., Sahlmann, J., & Bardalez Gagliuffi, D. 2016, ApJ, 827, 25, doi: 10.3847/0004-637X/827/1/25
  * Burgasser et al. (2007) Burgasser, A. J., Reid, I. N., Siegler, N., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 427.  https://arxiv.org/abs/astro-ph/0602122
  * Butler et al. (1996) Butler, R. P., Marcy, G. W., Williams, E., et al. 1996, PASP, 108, 500, doi: 10.1086/133755
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Close et al. (2012) Close, L. M., Puglisi, A., Males, J. R., et al. 2012, ApJ, 749, 180, doi: 10.1088/0004-637X/749/2/180
  * Cottaar & Hénault-Brunet (2014) Cottaar, M., & Hénault-Brunet, V. 2014, A&A, 562, A20, doi: 10.1051/0004-6361/201322445
  * Da Rio et al. (2014) Da Rio, N., Tan, J. C., & Jaehnig, K. 2014, ApJ, 795, 55, doi: 10.1088/0004-637X/795/1/55
  * Da Rio et al. (2016) Da Rio, N., Tan, J. C., Covey, K. R., et al. 2016, ApJ, 818, 59, doi: 10.3847/0004-637X/818/1/59
  * Da Rio et al. (2017) —. 2017, ApJ, 845, 105, doi: 10.3847/1538-4357/aa7a5b
  * D’Antona & Mazzitelli (1994) D’Antona, F., & Mazzitelli, I. 1994, ApJS, 90, 467, doi: 10.1086/191867
  * D’Orazi et al. (2009) D’Orazi, V., Randich, S., Flaccomio, E., et al. 2009, A&A, 501, 973, doi: 10.1051/0004-6361/200811241
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8, doi: 10.3847/0067-0049/222/1/8
  * Feiden (2016) Feiden, G. A. 2016, A&A, 593, A99, doi: 10.1051/0004-6361/201527613
  * Fűrész et al. (2008) Fűrész, G., Hartmann, L. W., Megeath, S. T., Szentgyorgyi, A. H., & Hamden, E. T. 2008, ApJ, 676, 1109, doi: 10.1086/525844
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306, doi: 10.1086/670067
  * Foster et al. (2015) Foster, J. B., Cottaar, M., Covey, K. R., et al. 2015, ApJ, 799, 136, doi: 10.1088/0004-637X/799/2/136
  * Fregeau et al. (2002) Fregeau, J. M., Joshi, K. J., Portegies Zwart, S. F., & Rasio, F. A. 2002, ApJ, 570, 171, doi: 10.1086/339576
  * Gaia Collaboration et al. (2016) Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1, doi: 10.1051/0004-6361/201629272
  * Gaia Collaboration et al. (2023) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2023, A&A, 674, A1, doi: 10.1051/0004-6361/202243940
  * Ginsburg et al. (2019) Ginsburg, A., Sipőcz, B. M., Brasseur, C. E., et al. 2019, AJ, 157, 98, doi: 10.3847/1538-3881/aafc33
  * Gutermuth et al. (2009) Gutermuth, R. A., Megeath, S. T., Myers, P. C., et al. 2009, ApJS, 184, 18, doi: 10.1088/0067-0049/184/1/18
  * Hillenbrand (1997) Hillenbrand, L. A. 1997, AJ, 113, 1733, doi: 10.1086/118389
  * Hillenbrand & Hartmann (1998) Hillenbrand, L. A., & Hartmann, L. W. 1998, ApJ, 492, 540, doi: 10.1086/305076
  * Hsu et al. (2021a) Hsu, C.-C., Theissen, Christopher A.and Burgasser, A. J., & Birky, J. L. 2021a, SMART: The Spectral Modeling Analysis and RV Tool, v1.0.0, Zenodo, doi: 10.5281/zenodo.4765258
  * Hsu et al. (2021b) Hsu, C.-C., Burgasser, A. J., Theissen, C. A., et al. 2021b, ApJS, 257, 45, doi: 10.3847/1538-4365/ac1c7d
  * Husser et al. (2013) Husser, T. O., Wende-von Berg, S., Dreizler, S., et al. 2013, A&A, 553, A6, doi: 10.1051/0004-6361/201219058
  * Kainulainen et al. (2017) Kainulainen, J., Stutz, A. M., Stanke, T., et al. 2017, A&A, 600, A141, doi: 10.1051/0004-6361/201628481
  * Karnath et al. (2019) Karnath, N., Prchlik, J. J., Gutermuth, R. A., et al. 2019, ApJ, 871, 46, doi: 10.3847/1538-4357/aaf4c1
  * Kim et al. (2019) Kim, D., Lu, J. R., Konopacky, Q., et al. 2019, AJ, 157, 109, doi: 10.3847/1538-3881/aafb09
  * Kounkel et al. (2016) Kounkel, M., Hartmann, L., Tobin, J. J., et al. 2016, ApJ, 821, 8, doi: 10.3847/0004-637X/821/1/8
  * Kounkel et al. (2022) Kounkel, M., Stassun, K. G., Covey, K., & Hartmann, L. 2022, MNRAS, 517, 161, doi: 10.1093/mnras/stac2695
  * Kounkel et al. (2018) Kounkel, M., Covey, K., Suárez, G., et al. 2018, AJ, 156, 84, doi: 10.3847/1538-3881/aad1f1
  * Kounkel et al. (2019) Kounkel, M., Covey, K., Moe, M., et al. 2019, AJ, 157, 196, doi: 10.3847/1538-3881/ab13b1
  * Kroupa (2000) Kroupa, P. 2000, New A, 4, 615, doi: 10.1016/S1384-1076(99)00048-2
  * Krumholz et al. (2014) Krumholz, M. R., Bate, M. R., Arce, H. G., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 243, doi: 10.2458/azu_uapress_9780816531240-ch011
  * Lacour et al. (2011) Lacour, S., Tuthill, P., Amico, P., et al. 2011, A&A, 532, A72, doi: 10.1051/0004-6361/201116712
  * Lada & Lada (2003) Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57, doi: 10.1146/annurev.astro.41.011802.094844
  * Majewski et al. (2017) Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94, doi: 10.3847/1538-3881/aa784d
  * Martin et al. (2018) Martin, E. C., Fitzgerald, M. P., McLean, I. S., et al. 2018, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10702, Ground-based and Airborne Instrumentation for Astronomy VII, ed. C. J. Evans, L. Simard, & H. Takami, 107020A, doi: 10.1117/12.2312266
  * Mathew & Federrath (2021) Mathew, S. S., & Federrath, C. 2021, MNRAS, 507, 2448, doi: 10.1093/mnras/stab2338
  * Matus Carrillo et al. (2023) Matus Carrillo, D. R., Fellhauer, M., Boekholt, T. C. N., Stutz, A., & Morales Inostroza, M. C. B. 2023, MNRAS, 522, 4238, doi: 10.1093/mnras/stad1256
  * McLean et al. (2000) McLean, I. S., Graham, J. R., Becklin, E. E., et al. 2000, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4008, Optical and IR Telescope Instrumentation and Detectors, ed. M. Iye & A. F. Moorwood, 1048–1055, doi: 10.1117/12.395422
  * McLean et al. (1998) McLean, I. S., Becklin, E. E., Bendiksen, O., et al. 1998, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 3354, Infrared Astronomical Instrumentation, ed. A. M. Fowler, 566–578, doi: 10.1117/12.317283
  * Moehler et al. (2014) Moehler, S., Modigliani, A., Freudling, W., et al. 2014, The Messenger, 158, 16
  * Morales-Calderón et al. (2012) Morales-Calderón, M., Stauffer, J. R., Stassun, K. G., et al. 2012, ApJ, 753, 149, doi: 10.1088/0004-637X/753/2/149
  * Murray & Correia (2010) Murray, C. D., & Correia, A. C. M. 2010, in Exoplanets, ed. S. Seager, 15–23, doi: 10.48550/arXiv.1009.1738
  * Murray & Dermott (1999) Murray, C. D., & Dermott, S. F. 1999, Solar system dynamics (Cambridge University Press)
  * Offner et al. (2023) Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 275, doi: 10.48550/arXiv.2203.10066
  * Palla & Stahler (1999) Palla, F., & Stahler, S. W. 1999, ApJ, 525, 772, doi: 10.1086/307928
  * Parker et al. (2016) Parker, R. J., Goodwin, S. P., Wright, N. J., Meyer, M. R., & Quanz, S. P. 2016, MNRAS, 459, L119, doi: 10.1093/mnrasl/slw061
  * Price-Whelan et al. (2017) Price-Whelan, A. M., Hogg, D. W., Foreman-Mackey, D., & Rix, H.-W. 2017, ApJ, 837, 20, doi: 10.3847/1538-4357/aa5e50
  * Raghavan et al. (2010) Raghavan, D., McAlister, H. A., Henry, T. J., et al. 2010, ApJS, 190, 1, doi: 10.1088/0067-0049/190/1/1
  * Reggiani et al. (2011) Reggiani, M., Robberto, M., Da Rio, N., et al. 2011, A&A, 534, A83, doi: 10.1051/0004-6361/201116946
  * Robberto et al. (2013) Robberto, M., Soderblom, D. R., Bergeron, E., et al. 2013, ApJS, 207, 10, doi: 10.1088/0067-0049/207/1/10
  * Scally et al. (2005) Scally, A., Clarke, C., & McCaughrean, M. J. 2005, MNRAS, 358, 742, doi: 10.1111/j.1365-2966.2004.08617.x
  * Strand (1958) Strand, K. A. 1958, ApJ, 128, 14, doi: 10.1086/146511
  * Stutz (2018) Stutz, A. M. 2018, MNRAS, 473, 4890, doi: 10.1093/mnras/stx2629
  * Stutz & Gould (2016) Stutz, A. M., & Gould, A. 2016, A&A, 590, A2, doi: 10.1051/0004-6361/201527979
  * Theissen et al. (2022) Theissen, C. A., Konopacky, Q. M., Lu, J. R., et al. 2022, ApJ, 926, 141, doi: 10.3847/1538-4357/ac3252
  * Tobin et al. (2009) Tobin, J. J., Hartmann, L., Furesz, G., Mateo, M., & Megeath, S. T. 2009, ApJ, 697, 1103, doi: 10.1088/0004-637X/697/2/1103
  * Tokovinin (2014) Tokovinin, A. 2014, AJ, 147, 87, doi: 10.1088/0004-6256/147/4/87
  * van Dam et al. (2006) van Dam, M. A., Bouchez, A. H., Le Mignant, D., et al. 2006, PASP, 118, 310, doi: 10.1086/499498
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261, doi: 10.1038/s41592-019-0686-2
  * Vitrichenko et al. (2006) Vitrichenko, É. A., Klochkova, V. G., & Tsymbal, V. V. 2006, Astrophysics, 49, 96, doi: 10.1007/s10511-006-0011-5
  * Weigelt et al. (1999) Weigelt, G., Balega, Y., Preibisch, T., et al. 1999, A&A, 347, L15, doi: 10.48550/arXiv.astro-ph/9906233
  * Winters et al. (2019) Winters, J. G., Henry, T. J., Jao, W.-C., et al. 2019, AJ, 157, 216, doi: 10.3847/1538-3881/ab05dc
  * Wizinowich et al. (2006) Wizinowich, P. L., Le Mignant, D., Bouchez, A. H., et al. 2006, PASP, 118, 297, doi: 10.1086/499290
  * York et al. (2000) York, D. G., Adelman, J., Anderson, John E., J., et al. 2000, AJ, 120, 1579, doi: 10.1086/301513