Stellar Proper Motions in the Orion Nebula Cluster

TL;DR

This study measures proper motions of stars in the Orion Nebula Cluster using HST and Keck data, revealing cluster dynamics, star escape candidates, and supporting virial equilibrium consistent with star formation theories.

Contribution

It provides the first detailed proper motion analysis of 701 stars in the ONC, confirming virial equilibrium and identifying high-velocity sources and their origins.

Findings

Over 10 escaping star candidates identified.

Proper motion dispersions consistent with virial equilibrium.

High-velocity sources BN, x, and n analyzed for origin.

Abstract

The Orion Nebula Cluster (ONC) is the nearest site of ongoing massive star formation, which allows us to study the kinematics and dynamics of the region in detail and constrain star formation theories. Using HST ACS/WFPC2/WFC3IR and Keck II NIRC2 data, we have measured the proper motions of 701 stars within an $\sim6'\times6'$ field of view around the center of the ONC. We have found more than 10 escaping star candidates, concentrated predominantly at the core of the cluster. The proper motions of the bound stars are consistent with a normal distribution, albeit elongated North-South along the Orion filament, with proper motion dispersions of $(\sigma_{\mu,\alpha^*}, \sigma_{\mu,\delta}) = (0.83\pm0.02,\,1.12\pm0.03)$ mas yr$^{-1}$ or intrinsic velocity dispersions of $(\sigma_{v,\alpha^*}, \sigma_{v,\delta}) = (1.57\pm0.04,\,2.12\pm0.06)$ km s$^{-1}$ assuming a distance of 400 pc to the ONC. The cluster shows no evidence for tangential-to-radial anisotropy. Our velocity dispersion profile agrees with the prediction from the observed stellar + gas density profile from Da Rio et al. (2014), indicating that the ONC is in virial equilibrium. This finding suggests that the cluster was formed with a low star formation efficiency per dynamical timescale based on comparisons with current star formation theories. Our survey also recovered high-velocity IR sources BN, x, and n in the BN/KL region. The estimated location of the first two sources $\sim500$ years ago agrees with that of the radio source I, consistent with their proposed common origin from a multi-stellar disintegration. However, source n appears to have a small proper motion and is unlikely to have been involved in the event.

Tables

Table 1. Table 1: HST observation

Epoch	Date	${\alpha }_{J2000}$	${\delta }_{J2000}$	Instrument	Filter	Exp Time	Proposal ID	PI Name
	(YYYY-mm-dd)	(hms)	(dms)		(nm)	(s)
1	1995-12-15	5:35:15.45	-5:24:06.65	WFPC2	F547M	200.0	6056	Rubin
	…	…	…	…	…	…	…	…
	1995-10-03	5:35:13.79	-5:21:47.13	WFPC2	F791W	100.0	5976	O’Dell
	…	…	…	…	…	…	…	…
2	1998-11-02	5:35:00.46	-5:24:40.00	WFPC2	F547M	500.0	6666	Stauffer
	…	…	…	…	…	…	…	…
	1998-11-02	5:35:00.46	-5:24:40.00	WFPC2	F791W	300.0	6666	Stauffer
	…	…	…	…	…	…	…	…
3	2000-09-13	5:35:13.77	-5:21:47.14	WFPC2	F547M	50.0	8121	O’Dell
	…	…	…	…	…	…	…	…
4	2001-03-13	5:35:17.00	-5:23:27.00	WFPC2	F439W	160.0	8894	Beckwith
	…	…	…	…	…	…	…	…
5	2004-10-12	5:35:18.43	-5:22:12.62	ACS	F435W	420.0	10246	Robberto
	…	…	…	…	…	…	…	…
	2004-10-12	5:35:18.43	-5:22:12.62	ACS	F555W	385.0	10246	Robberto
	…	…	…	…	…	…	…	…
	2004-10-12	5:35:18.43	-5:22:12.62	ACS	F775W	385.0	10246	Robberto
	…	…	…	…	…	…	…	…
6	2005-04-05	5:34:56.37	-5:23:19.89	ACS	F435W	420.0	10246	Robberto
	…	…	…	…	…	…	…	…
	2005-04-05	5:34:56.37	-5:23:19.89	ACS	F555W	385.0	10246	Robberto
	…	…	…	…	…	…	…	…
	2005-04-05	5:34:56.37	-5:23:19.89	ACS	F775W	385.0	10246	Robberto
	…	…	…	…	…	…	…	…
7	2007-09-12	5:35:12.05	-5:23:27.00	WFPC2	F791W	40.0	11038	Biretta
	…	…	…	…	…	…	…	…
	2007-09-12	5:35:12.05	-5:23:27.00	WFPC2	F547M	40.0	11038	Biretta
	…	…	…	…	…	…	…	…
8	2015-02-25	5:34:56.37	-5:23:19.89	ACS	F775W	340.0	13826	Robberto
	…	…	…	…	…	…	…	…
9	2015-03-11	5:34:47.07	-5:17:29.31	WFC3	F130N	302.9	13826	Robberto
	…	…	…	…	…	…	…	…
	2015-03-11	5:34:47.07	-5:17:29.31	WFC3	F139M	302.9	13826	Robberto
	…	…	…	…	…	…	…	…
10	2015-10-23	5:35:18.43	-5:22:12.62	ACS	F775W	340.0	13826	Robberto
	…	…	…	…	…	…	…	…
11	2015-10-24	5:35:09.34	-5:29:55.00	WFC3	F130N	302.9	13826	Robberto
	…	…	…	…	…	…	…	…
	2015-10-24	5:35:35.72	-5:31:04.42	WFC3	F139M	302.9	13826	Robberto
	…	…	…	…	…	…	…	…

Table 2. Table 2: Keck AO NIRC2 Wide-field Camera Observation

Epoch	Date	Filter	Exp Time	Nindi	Nstack	Total int. time	FWHM	Strehl
	(YYYY-mm-dd)		(s)			(s)	(mas)
12	2010-10-31	He I B	27.15	53	9	1439	162.32	0.136
13	2014-12-11	He I B	27.15	133	9	36111	132.40	0.191

Table 3. Table 3: PM measurements

ID	${\alpha }_{J2000}$aaEpoch 2015.5.	${\delta }_{J2000}$aaEpoch 2015.5.	${\mu }_{{\alpha }^{\ast }}$	${ϵ}_{{\mu }_{{\alpha }^{\ast }}}$	${\mu }_{\delta }$	${ϵ}_{{\mu }_{\delta }}$	Ndet	$\Delta$t	F139	NotebbIn this column, 1 = proplyd morphology; 2 = Herbig-Haro object (Ricci et al., 2008); 3 = double star, previously known in the literature (Hillenbrand, 1997; Robberto et al., 2013); 4 = new double star, previously reported as single in the literature; 5 = double-star candidate, previously reported as single in the literature; 6 = double-star candidate, previously unreported.
	deg	deg	mas yr-1	mas yr-1	mas yr-1	mas yr-1		yrs	mag
1	83.85346518	-5.36602620	0.34	0.09	-0.71	0.05	4	16.9	14.57
2	83.83552846	-5.34779414	-0.20	0.16	0.68	0.55	3	16.3	12.12
3	83.77062885	-5.35257430	-0.43	0.06	2.12	0.18	3	11.0	13.25
4	83.77579375	-5.33816393	-0.41	0.28	-0.75	0.30	3	11.0	14.52
5	83.82285598	-5.38091090	1.38	0.16	1.64	0.05	5	16.4	13.35
6	83.79410989	-5.37909779	0.14	0.46	-1.53	0.06	3	16.3	12.11
7	83.81797732	-5.34033710	-1.00	0.13	0.01	0.45	3	16.3	12.93
8	83.85596135	-5.39258964	-0.79	0.11	-0.44	0.11	4	16.3	12.81
9	83.84994021	-5.41941865	0.12	0.81	0.98	0.15	4	14.6	13.04
10	83.79543932	-5.37954202	-0.04	0.27	1.00	0.40	5	19.4	13.28
…	…	…	…	…	…	…	…	…	…	…

Table 4. Table 4: PMs as a function of magnitude

F139M	${\sigma }_{\mu ,{\alpha }^{\ast }}$	${\sigma }_{\mu ,\delta }$	${\sigma }_{\mu ,r}$	${\sigma }_{\mu ,t}$	${\sigma }_{\mu ,1D}$	N
mag	mas yr-1	mas yr-1	mas yr-1	mas yr-1	mas yr-1
$\mathrm{\hspace{0.17em}\hspace{0.17em}9.50}-12.18$	$0.86\pm 0.07$	$1.08\pm 0.08$	$0.91\pm 0.07$	$1.02\pm 0.08$	$0.98\pm 0.05$	110
$12.18-13.00$	$0.87\pm 0.07$	$1.18\pm 0.09$	$1.09\pm 0.08$	$1.00\pm 0.07$	$1.04\pm 0.06$	110
$13.00-13.56$	$0.84\pm 0.07$	$1.11\pm 0.09$	$0.92\pm 0.07$	$1.05\pm 0.08$	$0.98\pm 0.06$	110
$13.56-14.45$	$0.85\pm 0.07$	$1.19\pm 0.09$	$1.01\pm 0.08$	$1.07\pm 0.08$	$1.03\pm 0.06$	110
$14.45-16.09$	$0.83\pm 0.06$	$1.12\pm 0.08$	$1.06\pm 0.08$	$0.94\pm 0.07$	$0.99\pm 0.05$	110
$16.09-18.38$	$0.67\pm 0.05$	$1.03\pm 0.08$	$0.79\pm 0.06$	$0.94\pm 0.07$	$0.87\pm 0.05$	110

Table 5. Table 5: PM dispersions as a function of distance

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

Equations

$$L_i({\mu }_i|\bar{\mu },{\sigma }_{\mu })=G({\mu }_i|\bar{\mu },\sqrt{{\sigma }_{\mu }^2+{ϵ}_i^2}),$$

$$P(\bar{\mu },{\sigma }_{\mu }|D)\propto \prod _i{L}_i({\mu }_i|\bar{\mu },{\sigma }_{\mu })p(\bar{\mu },{\sigma }_{\mu }),$$

$$Q_{data,i}=\frac{{\mu }_i-\bar{\mu }}{\sqrt{{\sigma }_{\mu }^2+{ϵ}_i^2}}$$

$$Q_{theo,i}=\sqrt{2}{\mathrm{erf}}^{-1}(2({\mathrm{r}}_{\mathrm{i}}-0.5)/n-1),$$

$$⟨v_e^2{⟩}^{1/2}=2⟨v^2{⟩}^{1/2},$$

$${\bar{\mu }}_{{\alpha }^{\ast }}=-0.04\pm 0.03\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1},$$

$${\bar{\mu }}_{\delta }=-0.05\pm 0.05\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1},$$

$${\sigma }_{\mu ,{\alpha }^{\ast }}=0.83\pm 0.02\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1},$$

$${\sigma }_{\mu ,\delta }=1.12\pm 0.03\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}.$$

$${\bar{\mu }}_r=-0.04\pm 0.04\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1},$$

$${\bar{\mu }}_t=0.00\pm 0.04\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1},$$

$${\sigma }_{\mu ,r}=0.97\pm 0.03\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1},$$

$${\sigma }_{\mu ,t}=1.00\pm 0.03\mathrm{mas}\mathrm{y}{\mathrm{r}}^{-1}.$$

$$x(t)=x_{2014.9}+{\mu }_x(t-2014.9),$$

$$y(t)=y_{2014.9}+{\mu }_y(t-2014.9),$$

$$s(t)=\sqrt[]{x^2(t)+y^2(t)}$$

$$s_{min}=\frac{|x_{2014.9}{\mu }_x-y_{2014.9}{\mu }_y|}{\sqrt[]{{\mu }_x^2+{\mu }_y^2}}$$

$$t_{min}=-\frac{x_{2014.9}{\mu }_x+y_{2014.9}{\mu }_y}{{\mu }_x^2+{\mu }_y^2}$$

$$s_{min}=2\text{farcs}82\pm 1\text{farcs}32,$$

$$t_{min}=1535\pm 29\text{yr}.$$

Full text

proper motions — stars: kinematics and dynamics — stars: formation — open clusters and associations: individual (Orion Nebula Cluster)

††facilities: HST(ACS, WFPC2, WFC3), Keck:II (NIRC2), Gaia††software: Astropy, align, emcee, hst1pass, img2xym_WFC, img2xymrduv, Matplotlib, MVN, R, StarFinder

\noaffiliation

1 Introduction

The Orion Nebula Cluster (ONC) region provides an exquisite opportunity to probe the process of massive star and cluster formation in detail. The ONC is very massive, with stellar masses ranging between $0.1$ and $50\mathrm{M}_{\odot}$ (Hillenbrand, 1997). The mean age of the ONC is 2.2 Myr with a spread of a few Myr  (Reggiani et al., 2011), which is consistent with the star formation activity lasting between 1.5 and 3.5 Myr. The ONC’s close proximity ($\sim 400$ pc) and high galactic latitude ($b\sim 19^{\circ}$, or $\sim 135$ pc from the Galactic plane) allows us to study individual protostars and the entire cluster in detail. This combination is beneficial because the foreground has low extinction ($A_{V}$ = 1.5 mag; O’Dell & Yusef-Zadeh, 2000) and contains very few stars. Also, the Orion Molecular Cloud has a very large extinction up to $A_{V}$ = 50-100 mag (Hillenbrand & Carpenter, 2000; Scandariato et al., 2011), which reduces background confusion. Therefore, the stars observed in this region of the sky are mostly ONC members (Jones & Walker, 1988). The ONC allows us to probe the mechanisms that drive massive star and cluster formation, which remains a challenging problem in astrophysics.

Currently, two main theories attempt to explain massive stellar birth and they mainly differ in how and when the mass is gathered to form the star. The first model, called the Turbulent Fragmentation model, suggests that nearly the entire mass of individual protostars is gathered at a pre-stellar stage and that further fragmentation is halted due to external pressures from turbulence, radiation, and other forms of feedback  (McKee & Tan, 2002, 2003). The Competitive Accretion model, alternatively poses that mass is gathered during the star formation process itself, with all protostars starting with a low mass and accreting a significant amount of their final mass as they move through the molecular cloud (Bonnell et al., 2001a, b). One way to discern which model is more applicable is to study the dynamics of star-forming regions. While the Turbulent Fragmentation model requires the turbulence to remain virial and the star formation rates per dynamical timescale to be low, the competitive accretion model favors a rapid collapse of the gas clump and highly efficient star formation (Krumholz et al., 2011, 2012). Comparing the dynamical age of a star-forming cluster, such as the ONC, to the age spread of its stellar population may thus facilitate estimation of the star formation rates and distinguish the two models.

The dynamical properties of the stars can also have a significant impact on the star formation efficiency. Certain interactions could produce explosive outflows that provide feedback to the surrounding molecular cloud. The nature and frequency of these interactions inform our understanding of the role that feedback plays in halting star and cluster formation, expelling gas, and setting the overall star formation efficiency within a molecular cloud. Such an explosive event has been discovered in the ONC to the northwest of the well known Trapezium cluster (e.g. Zapata et al., 2005, 2006; Henney et al., 2007). This region hosts the Kleinmann-Low (KL) Nebula and contains a well studied radio and infrared source known as the Becklin-Neugebauer object (BN; Becklin & Neugebauer, 1967); thus the region is referred to as the BN/KL region. Based on analysis of the gas motions, the explosion is highly energetic ($2-6\times 10^{47}$ erg) and expelled over a very wide angle (Kwan & Scoville, 1976; Gómez et al., 2008; Bally et al., 2011), traced by molecules with a broad range of velocities ($>$100 km s*-1*; Kwan & Scoville, 1976; Furuya & Shinnaga, 2009; Bally et al., 2017). The outflow is the brightest known source of shocked H2 emission, with over 100 molecular bow shocks (e.g. Allen & Burton, 1993; Stolovy et al., 1998; Colgan et al., 2007). Millimeter and submillimeter observations suggest that the event was likely driven by close dynamical interactions in a group of massive protostars, including BN and source I, that resulted in a violent ejection of material and the formation of a compact binary or a stellar merger (Bally & Zinnecker, 2005; Bally et al., 2017).

There have been several previous studies of dynamical interactions and proper motions (PMs) within the Orion Nebula, both in the optical and radio. Originally, Parenago (1954) determined PMs for stars in the Orion Nebula over a field of $\sim$ 9 square degrees. Later, a 77 year baseline survey was done by  van Altena et al. (1988) for 73 stars in the Orion Nebula. Jones & Walker (1988) then carried out a survey using deep red-optical plates taken over 23 years on the Lick Shane reflector, which included over 1000 stars within 15$\arcmin$ of the ONC. In the radio, Gómez et al. (2005) measured the PMs of 35 sources in the Orion Nebula using the Very Large Array, with additional measurements presented in Gómez et al. (2008) and Dzib et al. (2017). Most recently, Kuhn et al. (2019) estimated the velocity dispersion of the ONC using the PMs of 50 sources in Gaia Data Release 2 (DR2; Gaia Collaboration et al., 2018) within $\sim 10\farcm 0$ of the center of the cluster. The ONC has proven a challenging environment for measuring PMs, particularly in the very center. These previous studies are limited either by their lack of precision or small sample size.

Fortunately, we now have access to a long baseline of data on the ONC from the Hubble Space Telescope and high-resolution near-IR data from the Keck II telescope focusing on the BN/KL region. Using these data, we have increased the precision of PMs, which has allowed us to further learn about the kinematics in this nearest massive star forming area.

We present the observations and data used to construct a new PM catalog for the ONC in §2. The analysis process for extracting PMs for each star is detailed in §3. The results are given in §4 followed by a discussion of how these results compare to previous studies in §5. Also in §5 we briefly discuss the interaction of sources near the BN/KL region.

2 Observations and Data

We measured stellar PMs near the center of the Trapezium and BN/KL region using high-resolution optical and infrared (IR) images spanning $\sim$20 years. Our final PM catalog covers $\sim 6\times 6$ sq arcmin around the Trapezium. The images were obtained with different instruments on board the Hubble Space Telescope (HST) including the Advanced Camera for Surveys with the Wide Field Channel (ACS/WFC), the Wide Field Camera 3 IR detector (WFC3/IR), and the Wide Field and Planetary Camera 2 (WFPC2) as well as the Near-Infrared Camera 2 (NIRC2) of the W. M. Keck II 10-m telescope.

2.1 HST

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

The HST images were obtained from several different cameras. The ACS/WFC consists of two 2048$\times$4096 pixel CCD detectors. The plate scale is 50 mas per pixel, which corresponds to a $202\arcsec\times 200\arcsec$ field of view. The WFPC2 uses four 800$\times$800 pixel CCDs where three of them cover a 150″$\times$150″region (WF) and have a pixel scale of 100 mas pixel*-1*. The fourth CCD (PC) images a $34\arcsec\times 34\arcsec$ square field with a spatial scale of 56 mas pixel*-1*. The WFC3 IR channel uses a single 1024$\times$1024 pixel CCD detector with a plate scale of 130 mas pixel*-1*, corresponding to a $136\arcsec\times 123\arcsec$ field of view.

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

2.2 Keck AO

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

The two epochs of Keck AO observations covered nearly the same sky area with the wide-field cameras on NIRC2, which has a pixel scale of 39.686 mas pixel*-1* (Yelda et al., 2010) and field of view of 40″$\times$40″, in the same passband, He I b ($\lambda_{0}=2.06\mu$m, $\Delta\lambda=0.03\mu$m). The narrow-band filter allows us to avoid the saturation of bright sources such as BN. The images were mosaicked around the BN/KL region for a total areal coverage of 1.4 square arcminutes. Sky frames were taken interspersed with science observations in a dark region $\sim 15\arcdeg$ to the east. Sky observations were timed in such a way that the field rotator mirror angle was identical to that of the science exposures, which is necessary to accurately subtract thermal emission from the field rotator mirror in this band (Stolte et al., 2008).

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

3 Analysis

3.1 Astrometry

3.1.1 HST

For ACS/WFC and WFC3/IR data, we used pipeline-calibrated images with a suffix _flt, which were dark and bias corrected and have been flat-fielded. All images were downloaded between February and June 2018 from the Mikulski Archive for Space Telescopes (MAST)111http://archive.stsci.edu. To measure stellar positions and fluxes in each exposure, we adopted the FORTRAN code hst1pass222http://www.stsci.edu/$\sim$jayander/CODE/, an advanced version of img2xym_WFC software package for HST (Anderson & King, 2006). The hst1pass code runs a single pass of source finding and point-source function (PSF) fitting for each exposure and corrects the positions of stars using the geometric-distortion correction of Anderson & King (2006) for ACS/WFC and the WFC3/IR correction developed by J. Anderson.333http://www.stsci.edu/$\sim$jayander/STDGDCs/ For WFPC2 data, we used calibrated images with a suffix _c0f and analyzed with the FORTRAN code img2xymrduv (Anderson & King, 2003). This code is implemented similarly to hst1pass and corrects the positions of stars from the WFPC2 data based on the distortion correction of Anderson & King (2003).

Outputs from hst1pass and img2xymrduv include the distortion-corrected positions of stars, their R.A. and decl. based on the WCS information in the images’ header, instrumental magnitudes, and the quality (or q) of the detections. Sources with q close to 0 appear very stellar, while those with large q values are mostly cosmic ray impacts or artifacts of diffraction spikes. For our analysis, we apply a quality cut with the threshold of $0<q\leq 0.5$ to exclude such false positives and saturated sources. We also set the minimum flux limits to be 1300 electrons for the narrow filter F130N and 500 electrons for other medium/wide filters, high enough to distinguish between the detections of stars and background noise.

3.1.2 Keck AO

The Keck AO NIRC2 data were reduced through a standard pipeline originally developed for analysis of Galactic center images (Stolte et al., 2008; Lu et al., 2009). This process includes dark and flat-field correction, sky subtraction, masking of bad pixels and cosmic rays, and the application of the distortion solution, provided by H. Fu444http://homepage.physics.uiowa.edu/$\sim$haifu/idl/nirc2wide/. The images were then registered and drizzled using the IRAF/PyRAF modules xregister and drizzle. The images were all stacked into one final average image for each pointing. Additionally, each final image had three associated subimages that combined one third of the data that were used to estimate astrometric and photometric uncertainties.

We used the IDL package StarFinder (Diolaiti et al., 2000) on each of the final averaged images and each of the subimages for each pointing with the wide camera to determine precise pixel positions for stars within the field. StarFinder extracts a PSF from the image over several iterations from a set of stars, for which we selected a total of five to seven stars in a range of magnitudes $8<K<10$.

The star catalogs from each pointing were then matched with the corresponding catalogs from their subimages with the program align (Ghez et al., 2008), which crossmatches sources and solves for astrometric transformations. The final star catalogs only include stars that appeared in each of the subimages as well as the averaged images. Finally, false or low-quality detections were rejected from the final star catalogs based on the PSF correlation coefficients of individual stars provided by StarFinder. We applied a threshold of $corr\geq 0.8$ for the quality cut, which is strict enough to reject unlikely detections while retaining faint stars (Diolaiti et al., 2000).

3.2 Relative PMs

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

We constructed a new reference frame, building upon Gaia DR2 using our star catalogs from ACS/WFC F775W and WFC3 IR F139M in Epochs #8-11 as follows. First, we converted the celestial coordinates of stars in Gaia DR 2 into right-handed Cartesian coordinates x and y parallel to the R.A. and decl. directions (i.e. x=$\Delta\alpha\cos\delta$, y=$\Delta\delta$)555Higher-order projection such as Equation (1) in van de Ven et al. (2006) changes our PM measurements typically 1 $\mu$as yr*-1* or less., respectively. In order to compromise between the number of reference sources and the accuracy of the astrometry, we adopted only Gaia stars with measurement errors smaller than 60 mas. We cross-matched the stars in the HST F775W and F139M catalogs with those of the Gaia DR2 catalog within a radius of 60 mas and applied 2D second-order polynomials to transform the positions of stars in the input star catalogs to those of the reference stars. Although the minimum number of data points for a second-order polynomial fit is 6, we excluded frames that had less than 9 matches to ensure a good fit. The median position was calculated for each star from both Gaia DR2 and the transformed catalogs. During each polynomial fitting, we rejected 3-sigma outliers and repeated the same process until convergence into a final solution.

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

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

3.3 PM dispersion calculation

From the measured PM measurements, we derived the internal PM dispersion of the ONC using Bayes’s theorem and a multivariate normal distribution model. Assuming that each PM measurement is drawn from a Gaussian distribution with mean $\bar{\mu}$ and an intrinsic dispersion $\sigma_{\mu}$, the likelihood for the ith PM measurement $\mu_{i}\pm\epsilon_{i}$ is then given by

[TABLE]

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

[TABLE]

where $p(\bar{\mu},\sigma_{\mu})\equiv p(\bar{\mu})p(\sigma_{\mu})$ is the prior for the mean and the standard deviation. We adopt a flat prior for the mean and a “non-informative” Jeffreys prior for the standard deviation - i.e., $p(\sigma_{\mu})\propto\sigma_{\mu}^{-1}$ (see §VII in Jaynes, 1968, for justification)666Flat prior can be strongly “informative” for a scale parameter and bias the posterior probability distribution (see, e.g., §4.1 of Eriksen et al., 2008). Nevertheless, using a flat prior instead of the Jeffreys prior $\sigma_{\mu}^{-1}$ increases our estimates of PM dispersions in this paper only by $\sim 0.01$ mas yr*-1* or less except for the sample in Figure 3, by only $\sim 0.05$ mas yr*-1*, due to its small sample size. The choice of prior thus changes none of our concolusion in this paper.. As each star has PM measurements along R.A. ($\alpha$) and decl.($\delta$), we maximized the product of the posteriori for the two axes, i.e. $P_{\alpha,\delta}=P_{\alpha}P_{\delta}$ utilizing the Markov chain Monte Carlo (MCMC) Ensemble sampler emcee (Foreman-Mackey et al., 2013).

4 Results

We present the PMs for 701 stars centered on the Trapezium and the BN/KL regions (Figure 2), adding $\sim 500$ sources with precise PM measurements as compared to Gaia DR2 over the same region. Our catalog has a temporal baseline of $\sim 20$ years and extends the wavelength coverage to the near-IR. Throughout the remainder of the paper, we use the notations $\mu_{\alpha^{*}}$ and $\mu_{\delta}$ to denote projected PMs where $\mu_{\alpha^{*}}\equiv\mu_{\alpha}\cos\delta$.

4.1 Consistency with Gaia DR2 PMs

We compare our PM measurements with those in Gaia DR2. Due to the strong nebulosity, the astrometric measurements from Gaia DR2 have overall lower quality for stars around the ONC compared with other nebula-free regions. As a result, we adopt a generous quality cut for Gaia DR2 stars in order to compromise between the astrometric quality and the size of the comparison sample777The Gaia DR2 sources with astrometric_excess_noise=0, adopted by Kuhn et al. (2019) for accurate estimates of measurement uncertainty, are generally too bright and saturated in our HST images (see Figure 1).: astrometric_gof_al$<$16 and photometric_mean_g_mag$<$16. With this condition, we found fifteen matches between Gaia DR2 and our PM catalog. The outer panel of Figure 3 shows the difference in PMs along the R. A. and decl. axes, where the data points are concentrated around (0,0) within 1 mas yr*-1*. Figure 4 verifies that the PM vectors tend to point in similar directions with similar magnitudes.

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

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

4.2 High-velocity Stars

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

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

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

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

[TABLE]

[TABLE]

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

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

[TABLE]

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

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

4.3 Internal PM dispersions

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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Our PM dispersions agree with those found by Jones & Walker (1988); ($\sigma_{\mu,\alpha^{*}},\sigma_{\mu,\delta}$)=($0.91\pm 0.06,1.18\pm 0.05$), ($\sigma_{\mu,r},\sigma_{\mu,t}$)=($1.06\pm 0.05,1.04\pm 0.05$) mas yr*-1*.

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

5 Discussion

5.1 Normality and isotropy of PM distribution

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

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

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

5.2 Dynamical equilibrium

It has been argued in some previous studies that the ONC is likely to be supervirial (e.g., Jones & Walker, 1988; Da Rio et al., 2014). The dynamical mass of the ONC inferred from the previous kinematic analysis done by Jones & Walker (1988) is nearly twice the total stellar and gas mass. Alternatively, virial equilibrium requires a one-dimensional mean velocity dispersion $\sigma_{v}\simeq 1.7\pm 0.3$ km s*-1* given the volume density of the stellar and gas contents in the ONC (Da Rio et al., 2014), which is only $\sim 75\%$ of the velocity dispersion of $2.34\pm 0.06$ km s*-1* from Jones & Walker (1988), leading to the conclusion that the ONC is likely to be slightly supervirial, with a virial ratio of $q\simeq 0.9\pm 0.3$.

This past result is, however, partially attributable to the previous estimate of the distance of the ONC adopted for the derivation of velocity dispersion in Jones & Walker (1988), $\sim 470$ pc.

The estimate of the distance has subsequently decreased in later studies to $\sim 400$ pc (e.g., Menten et al., 2007; Jeffries, 2007; Kounkel et al., 2017, 2018; Kuhn et al., 2019; Großschedl et al., 2018). At a distance of $d=414\pm 7$ pc from Menten et al. (2007), Jones & Walker (1988) would have obtained a smaller value for their velocity dispersion $\sigma_{v}\simeq 2.1\pm 0.1$ km s*-1*, which would give a virial ratio as low as $q\simeq 0.7\pm 0.3$.

Figure 8 compares our measured velocity dispersions in Table 5 to the predicted velocity dispersions required for virial equilibrium based on the total mass profile from Da Rio et al. (2014). We note that we adopted $414\pm 7$ pc from Menten et al. (2007) for the distance of the ONC as Da Rio et al. (2014) did, for consistency purposes. Overall, our velocity dispersion profiles tend to decrease with radius following the prediction based on the observed mass profile, a power-law profile slightly steeper than a singular isothermal sphere (Da Rio et al., 2014). When we include the escaping star candidates, our measured velocity dispersion appears to be more than $1\sigma$ larger than predicted at the very central region. However, when neglecting the escaping stars in ESC group 1 or 1 + 2, our measurements are in good agreement with the predicted values within $1\sigma$ uncertainty, suggesting that the bulk of the cluster is virialized. Even if the escaping star candidates are included, the measured velocity dispersions are still well below the boundedness limit ($\sigma_{bound}=\sqrt{2}\sigma_{vir}$). There is no indication of global expansion in the ONC from our analysis; the mean PM along the radial axis is small and consistent with zero within the $\sigma$ uncertainty, as shown in §4.3.

Kuhn et al. (2019) has reported evidence of mild expansion in the ONC based on PM measurements of Gaia DR2 sources, finding the median outward velocity $\widetilde{v_{out}}=0.42\pm 0.20$ km*-1* based on the uncertainty-weighted median with bootstrap resampling. In Figure 5 of their paper, the weighted kernel-density estimate (KDE) plot of $v_{out}$ exhibits two peaks, one at $\sim-0.4$ km s*-1* and the other at $\sim 1.0$ km s*-1*, which implies a significant concentration of data points or weights at two different places. To explore this, we divided the Gaia DR2 sample of Kuhn et al. (2019), including 378 stars, into two subgroups in terms of weights: (a) 42 stars with small errors (i.e. large weights) $\epsilon_{v,out}<0.15$ km s*-1* and (b) 336 stars with larger errors $\epsilon_{v,out}>0.15$ km s*-1*. The sum of the weights of the group (a) reaches $\sim 64$% of that of the group (b). For the groups (a) and (b), the median velocity is found to be $\widetilde{v_{out}}=0.94\pm 0.41$ km $s^{-1}$ and $0.09\pm 0.19$ km $s^{-1}$, respectively. We note that we converted PMs into $v_{out}$ and calculated the median velocity in the same manner as described in §3 and §4 of Kuhn et al. (2019), for consistency purposes. This result indicates that the $\sim 10$% of the whole sample (i.e. group (a)) produces the second peak in the KDE plot and biases the median velocity, ultimately leading to the conclusion that the ONC shows evidence for mild expansion. In fact, the stars in group (a) mostly fall into the magnitude range between phot_g_mean_mag $\sim 13$ and 15, which corresponds to the transition point in the Gaia error terms from the detector and calibration-dominated regime to photon noise-limited regimes. In this interval, the astrometric uncertainties of Gaia DR2 are known to be the most underestimated (see §4.6.4 and Figure 24 in Arenou et al., 2018). We refer to the IAU GA presentation slides by Lindegren888https://www.cosmos.esa.int/web/gaia/dr2-known-issues for more details. Between the group (a) and our catalog, we found one star matched, whose PM errors in the Gaia DR2 need to be increased by a factor of $\sim 3$ or greater for its outward velocities to be consistent within 1$\sigma$ uncertainty (see also §4.1). We therefore conclude that the net outward velocity claimed by Kuhn et al. (2019) is likely a result of underestimated uncertainties in the Gaia DR2 data.

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

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

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

5.3 Origin of high-velocity stars

Our PM catalog includes previously known fast-moving sources around the BN/KL complex, namely BN, source I, and source n, although in our analysis source n exhibits a rather small PM in the rest frame of the ONC. These optically invisible objects are all detected in epoch 2010 and 2014 NIRC2 He I b ($\sim$ K-band) images999The astrometric measurements of the highly embedded object BN on the HST WFC3/IR F130N and F139M ($\sim$ J and H-bands) images were excluded from this analysis as it appears highly elongated and asymmetric on those images.. Before the recent discovery of source I, two possible scenarios were proposed to explain their origin:

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

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

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

6 Conclusions

Using HST ACS/WFC3IR data that span $\sim 20$ years and Keck II NIRC2 data obtained in 2010 and 2014, we obtained relative PMs of 701 stars within $\sim 3\farcm 0$ of the ONC. With the analysis of these PMs, we reach the following conclusions.

Excluding the kinematic outliers, the PMs of our sample are consistent with a multivariate normal distribution. With the refined sample, the calculated velocity dispersions are $(\sigma_{\mu,\alpha^{*}},\,\sigma_{\delta})=(0.83\pm 0.02,\,1.12\pm 0.03)$ mas yr*-1* and $(\sigma_{\mu,r},\,\sigma_{\mu,t})=(0.97\pm 0.03,\,1.00\pm 0.03)$ mas yr*-1*. These values agree with those in previous surveys (e.g., Jones & Walker, 1988); but have a factor of $\sim 2$ improved precision. 2. 2.

The PM distribution appears elongated north-south with an axis ratio of $\sigma_{\mu,\alpha^{*}}/\sigma_{\mu,\delta}=0.74\pm 0.03$, resembling the stellar distribution which is also elongated north-south with $b/a\approx 0.7$ in the central region. On the other hand, the radial and tangential PMs are consistent with tangential-to-radial isotropy, as indicated by the low deviation from isotropy ($\sigma_{\mu,t}/\sigma_{\mu,r}-1)=0.03\pm 0.04$. 3. 3.

Compared to the prediction from the total density profile, our velocity dispersion profile is in good agreement with a virialized state. This suggests that the star-forming region is dynamically evolved. 4. 4.

Our analysis recovered the fast-moving IR sources in the BN/KL region, including BN, I, and n. The PMs of BN and source I are consistent with previous measurements in the literature, whereas source n exhibits a relatively small PM, as previously seen at IR and millimeter wavelengths. The estimated locations of BN and source I when the closest separation took place agree with the initial position of the radio source I, implying the dynamical decay of a multiple system involving these three sources. 5. 5.

The majority of ESCs are concentrated around the core of the ONC, where their PM vectors mostly point outward. 6. 6.

Based on comparisons with current star formation theories, our result suggests that the ONC is forming stars with a low star formation efficiency per dynamical timescale.

Our analysis shows that high spatial resolution, near-IR coverage of the ONC is essential; HST WFC3/IR + Keck NIRC2 observations revealed a factor of $\sim 3$ more stars than the optical Gaia DR2 sources in the BN/KL region. In order to obtain a complete PM catalog of the most embedded and the lowest mass objects, additional observations over a sufficiently long time baseline with near-IR telescopes/instruments such as HST WFC3/IR and Keck NIRC2 and the next generation near/mid-IR telescopes such as JWST and WFIRST will be required. Alongside the PM analysis, ongoing spectroscopic surveys for stellar line-of-sight velocities around the ONC will enable determining 3-dimensional stellar velocity distribution (e.g., Theissen et al. in preparation).

Acknowledgments. The authors thank the anonymous referee for interesting suggestions and comments. We thank Andrea Bellini, Joshua Simon, and Anthony Brown for helpful discussions on Gaia data, Gaspard Duchene on visual binaries in the ONC, and Michael Kuhn for providing the list of Gaia DR2 source IDs. We also thank Steven Stahler and Lynne Hillenbrand for informative discussions and comments. J. R. L. and D. K. acknowledge support from NSF AST-1764218, HST-AR-13258, and HST-GO-13826. J. R. L., Q. K., D. K., and C. A. T. are supported by the NSF grant AST-1714816. This work is based on observations made with the NASA/ESA Hubble Space Telescope, obtained at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. 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. Some of 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. We acknowledge use of the Astropy (Astropy Collaboration et al., 2013), Matplotlib (Hunter, 2007), and R (R Core Team, 2018) software packages.

Appendix A Systematic Errors in Gaia DR2 PMs

In order to use the Gaia DR2 as a control sample (see §4.1), it is important to understand underlying systematic errors as well as random errors depending on the quality of individual measurements. Here, we demonstrate that the Gaia DR2 PMs may include systematic errors induced by the scanning law of the survey, by comparing to HST PMs in the globular cluster NGC 7078. Note that NGC 7078 provides a cleaner data set than the nebulous ONC. The top left and right panels of Figure 10 show the PM vector point diagrams of stars in common between Gaia DR2 and the HST PM catalog of Bellini et al. (2014). The distribution of Gaia DR2 PMs exhibits unexpected linear structures nearly perpendicular to one another. These structures still remain in the bottom left panel, where the HST PMs are subtracted from the Gaia DR2 PMs. The linear trends closely resemble the scanning footprints around the cluster, which can be traced by the positions of Gaia DR2 sources filtered based on the number of good observations, astrometric_n_good_obs_al, as shown in the bottom right panel. This suggests that the systematic errors reflect the scanning pattern and survey incompleteness.

We find that fainter stars with a smaller number of observations and visibility periods tend to be more strongly effected by systematics. Figure 11 compares obvious outliers along the structures, highlighted in red in the left bottom panel of Figure 10, to the rest in different parameter spaces. The outliers are mostly fainter than phot_g_mean_mag$\sim 16$ and have astrometric_n_good_obs_al values below 120 and visibility_periods_used values below 8. The black solid lines in the left panel mark the quality cut thresholds applied for our control sample in §4.1 based on magnitudes and goodness-of-fit (gof) statistic of individual sources. These criteria alone cannot completely rule out the possibility of systematic errors; there are still a few contaminants, circled in blue, inside the boundaries. The right panels show that increasing thresholds for astrometric_n_good_obs_al and visibility_periods_used can reduce the number of stars subjected to the systematic errors in a control sample. Filtering based on these two parameters may come at a cost, especially where the detection efficiency of Gaia is typically low like at the central region of globular clusters due to high crowding (Arenou et al., 2018).

Appendix B The minimum separation of BN and source x in the past

We have determined the closest approach distance between BN and source I based on our PM measurements in a similar approach to the method of Gómez et al. (2008). Assuming their PMs are linear, the relative positions of BN with respect to source I as a function of time $t$ are given by

[TABLE]

[TABLE]

where ($x_{2014.9}$, $y_{2014.9}$) and ($\mu_{x},\mu_{y}$) are relative positions for epoch 2014.9 and relative motions in the right-handed Cartesian coordinates (x=$\Delta\alpha\cos\delta$, y=$\Delta\delta$). The separation of the two sources as a function of time is then given by

[TABLE]

The minimum separation $s_{min}$ and the corresponding epoch $t_{min}$ are given by differentiating the separation and equaling to 0:

[TABLE]

[TABLE]

Given a relative position of ($x_{2014.9}$, $y_{2014.9}$)= ($-16706.2$, 14309.9) mas and PMs of ($\mu_{x},\mu_{y}$) = ($-34.0\pm 3.1$, $30.6\pm 2.4$) mas yr*-1* based on our measurements for BN and source I, we obtain

[TABLE]

[TABLE]