SIRIUS Project – V. Formation of off-centre ionized bubbles associated with Orion Nebula Cluster

ABSTRACT

Massive stars born in star clusters terminate star cluster formation by ionizing the surrounding gas. This process is considered to be prevalent in young star clusters containing massive stars. The Orion Nebula is an excellent example associated with a forming star cluster including several massive stars (the Orion Nebula Cluster, ONC) and a 2-pc-sized H ii region (ionized bubble) opening towards the observer; however, the other side is still covered with dense molecular gas. Recent astrometric data acquired by the Gaia satellite revealed the stellar kinematics in this region. By comparing these data with star cluster formation simulation results, we demonstrate that massive stars born in the ONC centre were ejected via three-body encounters. Further, orbit analysis indicates that θ² Ori A, the second massive star in this region, was ejected from the ONC centre towards the observer and is now returning to the cluster centre. Such ejected massive stars can form a hole in the dense molecular cloud and contribute to the formation of the 2-pc bubble. Our results demonstrate that the dynamics of massive stars are essential for the formation of star clusters and H ii regions that are not always centred by massive stars.

methods: numerical, stars: kinematics and dynamics, H ii regions, open clusters and associations: individual: Orion Nebula Cluster

1 INTRODUCTION

The Orion Nebula (M42) is a well-known nebula in the Orion molecular cloud and is accompanied by a prominent ionized bubble (H ii region, see Fig. 1a). The 2-pc bubble (Extended Orion Nebula) is open towards the observer and considered to be ionized by massive stars in the Orion Nebula Cluster (ONC). The ONC is a young (∼1 Myr), massive (⁠|${\sim}2000\, \mathrm{M}_{\odot }$|⁠), and the closest (388 ± 5 pc from us) open cluster located in the brightest region of the Orion Nebula (Hillenbrand & Hartmann 1998; Palla & Stahler 1999; Kounkel et al. 2017). The trapezium, a group of massive stars in the dense core of the ONC, illuminates the surrounding gas and forms the H ii region (Güdel et al. 2008; O’Dell et al. 2009). In particular, θ¹ Ori C, an O-type star with a mass of 45 M_⊙ (Simón-Díaz et al. 2006), is the most energetic star in the ONC and is, therefore, a candidate source of energy for the formation of the 2-pc bubble (Güdel et al. 2008; Pabst et al. 2019; Geen et al. 2021).

Figure 1.

Panel (a): distribution of the O-type (blue circles) and B-type (green squares) stars. The three O-stars are NU Ori (HD 37061), θ¹ Ori C, and θ² Ori A from top to bottom. The line segment corresponds to the tangential velocity of the stars relative to the ONC (red cross). The physical scales x and y in pc are computed assuming that the ONC is 400 pc away from the Sun. The diffuse component in the background image (red band in the Digitized Sky Survey 2) shows the distribution of gas around the ONC region. Panel (b): the gas surface density and star distribution in the simulation (Model A). White, cyan, and magenta dots show stars with masses of >2, >10, and >20 M_⊙, respectively. Panel (c): the gas temperature of the simulation Model A in the range −1 < z < 3 pc, where the cluster is located at z = 0. The cyan and red dots show stars with masses of >2 and >20 M_⊙, respectively. The grey arrows indicate the velocities of massive stars.

As seen in the ONC, massive stars are closely related to the formation of star clusters. It is considered that the formation of star clusters is terminated by the formation of massive stars within them (Krumholz, McKee & Bland-Hawthorn 2019, and references therein). In previous star cluster formation simulations, however, the kinematics of massive stars has not been fully included. The gravitational forces between the stars were softened using a softening length (e.g. Bonnell, Bate & Vine 2003; Fukushima & Yajima 2021). This treatment prevents close encounters of stars and ejections of stars due to three-body encounters. Dinnbier & Walch (2020) performed simulations of star clusters embedded in a gas cloud including accurate stellar dynamics without gravitational softening. They demonstrated that star clusters are ionized from the outskirts of the clusters due to the dynamical scattering of massive stars in the cluster centre. They also reported that this process results in a change in the ionization time-scale. Another important influence shown in their simulations is that the ionization bubbles are off-centred when they resolve individual stars and integrate their accurate motions. Thus, the dynamics of massive stars may play an important role in the formation of star clusters.

Such an off-centre bubble is also seen in the Orion Nebula. It is considered that the 2-pc ionized bubble associated with the ONC is formed by θ¹ Ori C, which is the most massive star in this region and is located at the centre of the ONC. The expansion of the ionized bubble is estimated as

$$\begin{eqnarray*} R_{\rm s}(t)=[125/(154\pi)]^{1/5}(L_{\rm w}/\rho)^{1/5} t^{3/5}, \end{eqnarray*}$$

(1)

where ρ is the initial gas density, L_w is the mechanical luminosity of the stellar wind, and t is time (Weaver et al. 1977). Assuming the mechanical luminosity of θ¹ Ori C as L_w = 7.5 × 10³⁵ erg s⁻¹ and the density of the molecular cloud as n ∼ 10³ cm⁻³, Pabst et al. (2019) estimated that the ionized bubble has expanded to 2.2 pc in 0.2 Myr. The 2-pc H ii region is open towards the observer, and an expanding shell with a speed of ∼10 km s⁻¹ is observed (Pabst et al. 2019). Thus, θ¹ Ori C is estimated to be energetic enough to form the 2-pc bubble.

Another O-type star observed in the 2-pc bubble is θ² Ori A, which is the second massive star located ∼0.3 pc from the ONC centre in the projection (O’Dell, Kollatschny & Ferland 2017). The mass is estimated to be 25–39 M_⊙ (Preibisch et al. 1999; Simón-Díaz et al. 2006), and some protoplanetary discs in this region are suggested to be ionized by θ² Ori A rather than θ¹ Ori C (O’Dell et al. 2017). Although its mechanical luminosity is an order of magnitude smaller than that of θ¹ Ori C, θ² Ori A could have contributed to the formation of the current structure of the Orion Nebula as seen in the simulation of Dinnbier & Walch (2020). Furthermore, the NU Ori, which is an O-star located at the centre of M43 bubble, may also be ejected from the ONC. The dynamical structures of the ONC must be confirmed using the astrometric data in this region.

Recent data releases of Gaia (Gaia Collaboration et al. 2021) have revealed the dynamical activities of the ONC. Several runaway (>30 km s⁻¹) and walkaway (5–30 km s⁻¹) stars that can be ejected from the ONC have been found (McBride & Kounkel 2019; Platais et al. 2020; Schoettler et al. 2020). In our previous paper (Fujii et al. 2021a, hereafter Paper I), we performed a simulation of the star cluster formation using our new code asura+bridge (SIRIUS project; Fujii et al. 2021a,b,c; Hirai, Fujii & Saitoh 2021), with which we can precisely integrate the motion of stars without gravitational softening. In Paper I, we performed a simulation of a star cluster, which resembles the ONC, and showed its star formation history via hierarchical mergers. We also showed that the numbers of runaway and walkaway stars found around the ONC can be explained using the dynamical ejection of stars in the ONC.

In this study, we perform a series of simulations for the formation of star clusters starting from a turbulent molecular cloud with asura+bridge. We present the dynamical structure of massive stars inside star clusters and compare it with that of the ONC. We also investigate the motions of massive stars in the Orion Nebula using Gaia data.

2 METHODS

2.1 Numerical simulation

We performed the simulations using a smoothed-particle hydrodynamics (SPH) code, asura+bridge (Fujii et al. 2021c), in which the orbits of stars are integrated with a high-order integrator as collisional systems. asura+bridge combines the hydro part with the stellar part using a Hamiltonian splitting scheme, Bridge (Fujii et al. 2007). In our simulation, all the star particles are integrated using a particle–particle particle–tree (P³T) scheme (Oshino, Funato & Makino 2011; Iwasawa, Portegies Zwart & Makino 2015) using a code for star cluster simulation, petar (Wang et al. 2020b; Fujii et al. 2021b). In petar, a slow-down algorithmic regularization scheme (Wang, Nitadori & Makino 2020a) is included for binaries, and therefore tight binaries can be integrated accurately. asura+bridge includes an H ii region model (Fujii et al. 2021c), in which the Strömgren radius around a massive star is calculated using the gas distribution around the star. The Strömgren radius is iteratively determined to satisfy the balance between the photon count emitted from the massive star and the amount of gas absorbing the photons within the radius. The thermal feedback is given to the gas particles within the radius to maintain the gas temperature of 10⁴ K. Mechanical feedback owing to the stellar wind, which is correlated with the luminosity (Renaud et al. 2013), is also given to the gas particles within the radius as radial velocities.

We adopted initial conditions, which are the same as that used in Paper I (Model A) and smaller cloud model (Model B). Herein, we briefly summarize our initial set-up of the molecular cloud. We also adopted a homogeneous spherical molecular gas cloud with a turbulent velocity field (Bonnell et al. 2003). For Model A, we set a mass of 2 × 10⁴ M_⊙ and a radius of 12 pc. The gas-particle mass is set as 0.01 M_⊙. Under this condition, the resulting initial gas density and free-fall time are 79.9 cm⁻³ and 4.87 Myr, respectively. For Model B, we adopt a mass of 5 × 10³ M_⊙ and a radius of 2.285 pc. The resulting density and free-fall time are 2.89 × 10³ cm⁻³ and 0.81 Myr, respectively. We adopted 0.1 M_⊙ for the gas-particle mass for Model B to reduce the calculation cost. The gas-mass resolution of 0.1 M_⊙ can result in an increase in the formation of stellar mass up to ∼10  per cent, but it does not affect the structures of star clusters (Fujii et al. 2021c).

A turbulent velocity field with power spectrum of ∝v⁻⁴ was given to the initial condition. We set the initial virial ratio (kinetic energy/potential energy) as 0.45 and 0.5 for Models A and B, respectively. The turbulent velocity is scaled to satisfy the given virial ratio. These models are similar to that of Bonnell et al. (2003). Because the result depends on the randomness of the initial turbulent velocity field, we performed multiple runs with different random seeds for the turbulence for each model. We performed three and five runs for Models A and B, respectively. We summarize these parameters in Table 1.

Table 1.

Models and parameters for the simulations of star cluster formation.

Name	Mg	mg	Rg	nini	tff, ini	αvir	ϵg	ϵs	nth	rmax	ΔtB	Nsin	Nrun
	(M⊙)	(M⊙)	(pc)	(cm−3)	(Myr)		(pc)	(pc)	(cm−3)	(pc)	(yr)
Model A	2 × 104	0.01	12	79.9	4.87	0.45	0.07	0.0	7.4 × 104	0.2	200	1	3
Model B	5 × 103	0.1	2.285	2.89 × 103	0.81	0.5	0.07	0.0	7.4 × 104	0.2	200	4	5

Note. From the left: model name, initial cloud mass (M_g), gas-particle mass (m_g), initial cloud radius (R_g), initial cloud density (n_ini), initial free-fall time (t_{ff, ini}), initial virial ratio (α_vir = |E_k|/|E_p|), softening length for gas (ϵ_g) and stars (ϵ_s), star formation threshold density (n_th), the maximum search radius (r_max), time-step for Bridge (Δt_B), the number of runs that resulted in the formation of a single star cluster (N_sin), and the number of runs (N_run).

For star formation, we used a probabilistic method often used for galaxy simulations; however, we modified it for star formation resolving individual stars (Hirai et al. 2021). In our scheme, stars are formed when gas particles satisfy the following conditions: (1) the gas density exceeds the threshold density (n_th), (2) the gas temperature is lower than the threshold temperature (T_th), and (3) the divergence of the velocity is less than zero. Here, we assume n_th = 7.4 × 10⁴ cm⁻³ and T_th = 20 K. We selected gas particles to star particles following Schmidt’s law (Schmidt 1959). We thereafter assigned stellar mass following an initial mass function that we assumed. We adopt the Kroupa mass function (Kroupa & Weidner 2003). As the gas-particle mass is always smaller than stellar masses, the new forming stars are created by assembling mass from a given search radius (r_max = 0.02 pc). Masses of gas particles within r_max are transported to newly formed stars. We limited the stellar mass to be half of the gas mass within r_max. If the chosen stellar mass exceeds this limit, we re-assigned the stellar mass until these conditions are satisfied.

Once star particles are created, we assigned the position and velocity of star particles to ensure momentum conservation. Once the mass of gas particles becomes 10 times smaller than its original mass, we merge the gas particle to the nearest gas particles to ensure that the number of gas particles is reasonably low. By using this model, we confirmed that it can reproduce the relationship between maximum mass and the enclosed stellar mass (Weidner, Kroupa & Pflamm-Altenburg 2013). This star formation procedure is implemented using a chemical evolution and star formation library, CELib (Saitoh 2017). These schemes are all tested in previous papers (Fujii et al. 2021b,c; Hirai et al. 2021).

The gravitational forces between gas particles and between gas and star particles are softened with a Plummer-type softening, and the softening length (ϵ_g) is set to be 0.07 pc. However, we did not use any softening for star particles. This treatment allows dynamical formation of binaries and scattering of stars owing to three-body encounters. Every 200 yr (Bridge time-step, Δt_B), stars are perturbed by the gravity of the gas. Moreover, gas particles are integrated using a tree algorithm with hierarchical time-steps.

We performed the simulation up to 10 Myr for Model A and ∼4 Myr for Model B, at which the molecular cloud is fully ionized.

2.2 OB stars near the ONC from Gaia catalogue

To compare the simulation with the ONC, we compiled a sample of OB stars near the ONC using Gaia Early Data Release 3 (EDR3) in two steps. First, we selected all the Gaia sources that satisfy the following properties: (i) Sky position is within 10° from the ONC; (ii) Gaia’s G-band magnitude is brighter than 12 mag; and (iii) Gaia’s parallax measurement ϖ ± σ_ϖ is consistent with a distance of 340–460 pc within 3σ. Namely, ϖ + 3σ_ϖ < 1/(340 pc) and ϖ − 3σ_ϖ > 1/(460 pc). Next, we cross-matched this sample against the Set of Identifications, Measurements and Bibliography for Astronomical Data (SIMBAD) catalogue and selected all the stars that were classified as O- or B-type stars in SIMBAD. This procedure resulted in 353 OB stars within 10° (corresponding to 70 pc at d = 390 pc) from the ONC, 34 of which were within 0|${_{.}^{\circ}}$|71 (5 pc) from the ONC.

The tangential velocity of OB stars relative to the ONC is given as (Δv_α, Δv_δ) = (4.74047 km s⁻¹) × (d/pc) × [(μ_α*, μ_δ) − (μ_α*, μ_δ)_ONC]/(mas yr⁻¹). Here, the heliocentric distance of each star is estimated to be d = 1/ϖ, where ϖ is the parallax determined by Gaia. Most OB stars in our sample had a good parallax measurement; hence, this simple point estimate of the distance (ignoring the uncertainty) is sufficient for our study.

Using the proper motion, we calculated the velocities relative to the ONC. The proper motion of the ONC’s centre of mass was assumed to be (μ_α*, μ_δ)_ONC = (1.1, 0.3) mas yr⁻¹ (Jerabkova et al. 2019). We summarized the proper motions and velocity relative to the ONC of the 34 OB stars in Table A1 in Appendix A.

Table 2.

Simulated star clusters.

Name	tONC, sim	tONC	Ms, 3pc	Mg, 3pc	vesc
	(Myr)	(Myr)	(M⊙)	(M⊙)	(km s−1)
Model A	7.7	3.2	5200	4540	5.29
Model B1	2.3	1.6	1590	1580	3.01
Model B2	2.45	1.75	1630	1200	2.85
Model B3	2.4	1.7	2000	1990	3.38
Model B4	2.4	1.7	1630	1640	3.07

Note. From the left: model name, time at which we compare the model with the ONC from the beginning of the simulation and the star formation (t_{ONC, sim} and t_ONC), stellar and gas masses within 3 pc from the cluster centre (M_{s, 3pc} and M_{g, 3pc}), and escape velocity at 3 pc (v_esc). The numbers in the model names indicate the runs with different random seeds for the initial turbulent velocity field.

Table 3.

Stellar parameters.

Name	Teff	log L/L⊙	M	|$v_{\inf }$|	|$\dot{M}$|	Lw
	(K)		(M⊙)	(km s−1)	(M⊙ yr−1)	(erg s−1)
θ1 Ori C	39 000	5.31	45	2500	4.0 × 10−7	7.9 × 1035
θ2 Ori A	35 000	4.93	39	2000	5.9 × 10−8	7.4 × 1034
NU Ori	31 000	4.42	16	1000	1.0 × 10−8	3.1 × 1033

Note. From the left, column 1 shows the name of the star, and columns 2, 3, and 4 represent the effective temperature, luminosity, and mass (Simón-Díaz et al. 2006). Column 5 shows the terminal velocity of the stellar wind (Prinja, Barlow & Howarth 1990; Nebot Gómez-Morán & Oskinova 2018). Column 6 shows the mass-loss rate obtained using equation (12) of Vink, de Koter & Lamers (2000) with the effective temperature (T_eff), luminosity (L), mass (M), and the ratio of terminal velocity (v_inf) and escape velocity (v_esc). We assumed that v_inf/v_esc = 2.6. Column 7 represents the mechanical luminosity of the wind, |$L_{\rm w}=0.5\dot{M}v_{\rm inf}^2$|⁠.

3 RESULTS

3.1 Star cluster formation in the simulation

In Fig. 2, we present the time evolution of gas surface density and temperature distribution with the distribution of stars. The star formation begins after about one initial free-fall time, which are at 4.87 and 0.81 Myr for Models A and B, respectively. Hereafter, we indicate the time from the beginning of the simulation as t_sim and the time from the first star formation (4.5 and 0.7  Myr for Models A and B, respectively) as t. The cluster evolves via mergers of stellar clumps and the outer region is ionized first owing to the dynamical ejection of massive stars.

Figure 2.

Snapshots of the simulation (Model A). Left-hand column: gas surface density. Dots indicate stars with >1 M_⊙. Cyan and red indicate stars with |$10\lt m\lt 20\, \mathrm{M}_{\odot }$| and >20 M_⊙, respectively. Right-hand column: gas temperature. Cyan dots indicate stars with >1 M_⊙. Red indicates stars with >20 M_⊙. Time indicates the time from the beginning of the simulation.

The total stellar mass evolution in the simulation is shown in Fig. 3. We assumed t = 3.2 Myr (t_sim = 7.7 Myr), at which the stellar and gas masses in the cluster are comparable, as the present day of the ONC for Model A. Meanwhile, a part of the gas in the cluster centre was ionized; however, some gas remained as a dense molecular cloud, and therefore star formation is still ongoing in the central region of the cluster. The age of stars in the ONC distributes within 10 Myr, but the majority were formed within the past 2 Myr and the star formation rate is accelerated (Palla & Stahler 1999). Our simulation is consistent with such observational features of the star formation in the ONC. The average stellar age at this time is 1.7 Myr, which is also consistent with the age of the ONC (1–3 Myr; e.g. Hillenbrand & Hartmann 1998; Da Rio et al. 2010).

Figure 3.

Stellar mass evolution in the simulations; Model A (top) and Model B (bottom). Star symbol indicates the time when we compared to the ONC (the time at which the gas and stellar masses in the cluster are comparable). The x-axis shows the time from the first star formation, which is 4.5 Myr (Model A) and 0.7 Myr (Model B) from the beginning of the simulation (top axis, t_sim).

The radial density distributions of gas and stars at t = 3.2 Myr are shown in Fig. 4. The central region within <0.3 pc has a density higher than 10⁵ cm⁻³, which is sufficient for star formation (Onishi et al. 2002). The cold gas density continues to increase at the centre of the cluster until 3.2 Myr. The central density of cold gas reaches 10⁷ cm⁻³, which is comparable to the observed molecular gas density in the ONC. The Atacama Large Millimeter/submillimeter Array N₂H⁺ (J= 1–0) and Atacama Pathfinder Experiment N₂H⁺ (J= 7–6) observations towards ONC revealed that the density of molecular clouds near the ONC centre is estimated to reach 10⁶ cm⁻³ and can partially exceed 10⁷ cm⁻³ (Hacar et al. 2020).

Figure 4.

Density profiles of stars and gas at t = 3.2 (t_sim = 7.7) Myr from the beginning of the star formation (simulation) for Model A. Black solid and dotted curves indicate the stellar and gas distributions, respectively. Blue and red dashed curves indicate cold gas with a temperature of <30 K and warm and hot gas with a temperature of >30 K.

After t ∼ 3.5 Myr, the star formation slows down (see Fig. 3), and at a later stage, the central cold gas density decreases. At the end of this simulation (t = 4.5 Myr), the amount of cold gas in the cluster centre is only ∼10 M_⊙, and it can no longer form new stars.

Once massive stars form in the simulation, the feedback ionizes and blows away the surrounding gas (see Fig. 2). However, it takes a few Myr for the massive stars to fully ionize the molecular cloud. In the panels (b) and (c) of Fig. 1, the snapshots at t = 3.2 Myr are shown. The total stellar mass of the simulated cluster at this time was 5200 M_⊙ within 3 pc from the cluster centre (see Table 2). In these snapshots, an ionized bubble similar to the Orion Nebula is observed, and the star formation is still ongoing in the dense molecular cloud close to the cluster centre. These structures are similar to those of the Orion Nebula (see Fig. 1a). This bubble is formed by a star ejected from the cluster centre due to strong dynamical interactions such as three-body encounters.

The mass of the star clusters that are formed and the number of star clusters change depending on the randomness of the initial turbulent velocity field. Although we performed three and five runs for Models A and B, respectively, one and four of them resulted in the formation of single star clusters, which is similar to the ONC. The other runs resulted in twin or multiple star clusters. We exclude these runs from the following analyses.

3.2 Velocity distribution of stars around the ONC

If the ONC has been formed via star formation in turbulent molecular clouds, stars formed in the central region of the cluster must be scattered from the cluster centre, and this mechanism must be confirmed by the velocity distribution of stars around the ONC. Therefore, we compared the velocity distribution of the stars within 5 pc from the ONC to that obtained from the simulation. The Gaia astrometry satellite has added kinematic data (positions and velocities) of stars in this region (see Fig. 1a), which provide insights into the dynamical evolution of this region (Gaia Collaboration et al. 2021).

Fig. 5 shows the cumulative velocity distribution of the massive stars (>2 M_⊙) within 5 pc from the clusters for the simulations and observations. We compare our simulations to the observations at the time when stellar and gas masses inside 3 pc from the cluster centre are comparable. The velocities of the stars with respect to the ONC are obtained from Gaia data (Gaia Collaboration et al. 2021). The total stellar mass of the ONC is estimated to be 1800–2700 M_⊙, and the total mass including the gas within ∼3 pc is estimated to be ∼5000 M_⊙ (Hillenbrand & Hartmann 1998; Tan, Krumholz & McKee 2006). Because the total masses of the clusters in the simulations are not exactly the same as that of the ONC, we scaled the number of stars and velocity relative to the cluster by the total mass and escape velocity of the cluster, respectively. The cluster and gas masses at the time we adopted are summarized in Table 2. While Model A was twice as big as the ONC, all of Model B were less massive.

Figure 5.

Velocity distribution of stars within 5 pc from the star cluster. The green curve represents the OB stars observed near the ONC (50th percentile). The grey curves indicate the 2.5, 16, 84, and 97.5 percentiles (see the text). The red dotted and magenta dashed curves indicate stars observed in our simulations between 0.8 and 5 pc from the centre of the formed cluster. Here, the velocities in the simulation are multiplied by |$\sqrt{2/3}$|⁠. The velocities of simulations are scaled by the escape velocities measured for both the observation and simulation. We used the total (gas + stellar) masses within 3 pc of |$5000\, \mathrm{M}_{\odot }$| (Hillenbrand & Hartmann 1998; Tan et al. 2006) for the observation. The cumulative numbers of stars are also scaled to the ONC mass (2500 M_⊙) for the simulations. The total mass and escape velocities of the simulations are summarized in Table 2. Thick and thin curves indicate Models A and B, respectively.

As the number of observed OB stars is only 34, we evaluate the Poisson noise in the cumulative tangential velocity distribution of the 34 OB stars near the ONC (see Table A1). First, we denoted the magnitude of the tangential velocity as v_{obs, n} for nth star in our sample (n = 1,…, 34), without taking into account the observational error. Then, we created 100 000 bootstrapped samples |$\lbrace v^{(j)}_i \mid i= 1,\dots ,34 \rbrace$|⁠, where j = 1,…, 100 000 denotes jth bootstrap sample. We note that ith star in jth bootstrap sample is randomly chosen from {v_{obs, n}}, and that we are allowed to choose the same star multiple times. For jth bootstrap sample, the cumulative distribution of |$\lbrace v^{(j)}_i \rbrace$| is expressed as CDF^(j)(v). For any given value of velocity v, we computed the 2.5, 16, 50, 84, and 97.5 percentile points of {CDF^(j)(v)∣j = 1,…, 100 000} using the 100 000 cumulative distribution functions. By varying v, we evaluated how these percentile points change as a function of v. We present the result in Fig. 5.

The velocity distributions of the massive stars were identical for both the simulations and observations. The number of stars drops around the escape velocity; however, the distribution continues beyond the escape velocity following a power law of −1.5. This power-law distribution is formed with high-velocity stars scattered inside the star cluster owing to three-body encounters (Perets & Šubr 2012). Thus, the velocity distribution of stars around the ONC suggests that scattering of massive stars occurred in the ONC.

Some stars were ejected from the cluster with a velocity higher than 30 km s⁻¹. They escape from the cluster and can be recognized as runaway stars. Slower stars (5–30 km s⁻¹) are called walkaway stars. Notably, some runaway and walkaway candidates originated from the ONC have been found (McBride & Kounkel 2019; Platais et al. 2020; Schoettler et al. 2020). We quantitatively discussed the runaway and walkaway stars around the ONC in Paper I.

3.3 Ionization due to scattered massive stars

The scattered massive stars contribute to the formation of the off-centre ionized bubbles. In Fig. 6, we show the distributions of massive stars and cold (<100 K) gas. The right-hand side of this panel is the z-direction in the right-hand panels of Fig. 1. As shown in this figure, a star cluster is formed in the centre of filament along the y-axis, and in the innermost region, massive stars form a 0.1-pc-scale H ii region.

Figure 6.

Edge-on density distribution of the molecular gas (T < 100 K) at the central region of the cluster in Model A for t = 3.0, 3.1, and 3.2 Myr from left to right. Panel (c) shows the same time for the snapshots as shown in Fig. 1. The z-axis is the direction of the observer. The star symbols are massive stars with >16 M_⊙ within x = ±0.3 pc. Yellow star is the most massive star in this region. Magenta and cyan stars are stars that are scattered from the cluster centre and then return.

Initially (panel a), the inner 0.1-pc bubble is covered with molecular gas. However, in panels (b) and (c), the small bubble is open towards the observer. This can be attributed to the ejection of a massive star from the central region. This structure is similar to that observed in the Orion Nebula (O’Dell et al. 2009; Abel, Ferland & O’Dell 2019; O’Dell, Abel & Ferland 2020). Behind the 2-pc bubble, a small 0.1-pc-scale bubble exists.

Another O-type star observed in the 2-pc bubble is θ² Ori A, which is the second massive star located ∼0.3 pc from the ONC centre in the projection (O’Dell et al. 2017). The mass is estimated to be 25–39  M_⊙ (Preibisch et al. 1999; Simón-Díaz et al. 2006), and some protoplanetary discs in this region are suggested to be ionized by θ² Ori A rather than θ¹ Ori C (O’Dell et al. 2017). We estimated the mechanical luminosity of θ² Ori A to be 7.4 × 10³⁴ erg s⁻¹ (see Table 3 and Appendix A). Although it is an order of magnitude smaller than that of θ¹ Ori C, the bubble expansion time to 0.1 pc driven by θ² Ori A is only 0.05 Myr for n = 10⁷ cm⁻³. This was obtained using equation (1). This expansion time is short enough for θ² Ori A to partially ionize the filament’s wall when it travels to the outer region due to the dynamical ejection.

We estimated the trajectory of θ² Ori A using its astrometric data. According to Gaia data, this star is located at a distance of |$336^{+26}_{-22}$| pc (1σ error) from the observer. This distance seems to be much closer to the observer than the 2-pc bubble, but we note that the parallax error of this star is ∼10  per cent, which is an order of magnitude larger than the other stars (see Table A1). Owing to 2σ error, this star is located at a distance of 295–391 pc. The velocity relative to the ONC is 3.3 km s⁻¹ towards the ONC (see Fig. 1). The proper motion suggests that θ² Ori A was ejected from the ONC centre and is currently returning to the ONC. θ² Ori A could break the molecular cloud towards the observer.

Assuming the mass distribution of the ONC, we calculate possible orbits of θ² Ori A. By fitting a double power-law function to the mass distribution in Fig. 4 and scaling the mass to the ONC mass, we obtained a mass distribution of |$M(\lt r)=25\,000\, (r/{\rm 1\,pc})^{1.3}\,\mathrm{M}_{\odot }$| for |$r\lt 0.175\, {\rm pc}$| and |$M(\lt r)=4000\, (r/{\rm 1\,pc})^{0.25}\,\mathrm{M}_{\odot }$| for |$r\gt 0.175\, {\rm pc}$|⁠. The enclosed mass at 3 pc was 5300 M_⊙, which is consistent with the observed ONC mass including gas (Hillenbrand & Hartmann 1998; Tan et al. 2006). We set the initial position as 0.01 pc and integrated the radial motions of stars by changing the initial radial velocity from 21.4 to 22.2 km s⁻¹. In this velocity range, the orbit satisfied 0.3 pc and −3 km s⁻¹ in projection within 0.1–0.4 Myr.

Some possible orbits (radial distances and velocities) are shown in Fig. 7. With the current distance of the θ² Ori A from the ONC centre of ∼0.3 pc and the velocity to the ONC of 3 km s⁻¹, we obtained an ejection time shorter than ∼0.4 Myr. This time-scale is consistent with the age of the 2-pc bubble, which was estimated from the shell velocity (∼0.2 Myr; Pabst et al. 2019, 2020).

Figure 7.

Possible orbital evolution of stars ejected from the ONC centre. Black and red curves indicate the distances from the ONC centre and the velocities relative to the ONC, respectively. The curves are for different initial velocities. Blue vertical lines indicate times (0.08, 0.17, and 0.38 Myr) at which the current projected distance (0.3 pc) and the current projected velocity (−3 km s⁻¹) of the θ² Ori A are satisfied. Stars indicate the positions and velocities at these times.

Once θ² Ori A was ejected from the ONC centre, it contributes to the ionization of the off-centre region of the ONC. Assuming the gas density of the outskirt of the cluster as 10³–10⁴ cm⁻³ (see Fig. 4), we obtained the bubble size of 1–2 pc for 0.3 Myr. In O’Dell et al. (2017), some protoplanetary discs are suggested to be ionized by θ² Ori A rather than θ¹ Ori C.

The structure of the ONC is shown in Fig. 8. θ² Ori A was ejected from the ONC <0.5 Myr ago and now is returning to the ONC centre. Through the hole opened by the ejection of θ² Ori A, θ¹ Ori C can irradiate the gas towards the observer (Güdel et al. 2008; Pabst et al. 2019). Both θ¹ Ori C and θ² Ori A contribute to the formation of the 2-pc bubble as is suggested from the observations of this region (O’Dell et al. 2017).

Figure 8.

Schematic of the Orion Nebula structure. The grey region indicates a dense Orion molecular cloud. The star symbols indicate the O-type stars, and orange symbols indicate other stars. The blue spheres indicate the H ii regions.

We also calculate the bubble size and the trajectory of another O-star in the Orion Nebula, NU Ori, which is located in the centre of M43 (see Fig. 1). The mass of this star is estimated to be 16 ± 3 M_⊙ (Gravity Collaboration et al. 2018), and the bubble radius around it is ∼0.3 pc. The bubble is still embedded in the dense molecular cloud. The velocity relative to the ONC is 2.8 km s⁻¹ in projection, and therefore we can estimate that this star has travelled from the ONC centre. From the velocity and distance, we estimated that this star was ejected 0.1–0.2 Myr ago. The distance to NU Ori (⁠|$415^{+11}_{-10}$| pc) suggests that this star appears to be behind the ONC, where the dense molecular gas still exists. Using equation (1), we estimated that the size of 0.3 pc at 0.1 Myr is achieved for the wind luminosity of NU Ori (L_w ∼ 10³³ erg s⁻¹) and the gas density of 10⁴ cm⁻³. These results are also consistent with observations.

4 SUMMARY

We performed N-body/SPH simulations of star cluster formation. Our simulation included the proper integration of stellar orbits, which revealed that star clusters eject massive stars from the cluster centre. Such ejected massive stars break the wall of the dense molecular cloud in the direction they escaped. The hole in the dense molecular cloud results in the formation of off-centre H ii regions such as the H ii region associated with the Orion Nebula. In contrast, the molecular gas is too dense to be completely ionized, and therefore the star formation is still ongoing.

We found that the distribution of velocities of OB stars around the ONC obtained from Gaia astrometric data is identical to the velocity distribution obtained from our simulations. This result implies that these OB stars around the ONC are formed in the ONC and scattered in the ONC centre. Some scattered massive stars did not escape from the cluster and fell back to the cluster centre. The astrometric data obtained from Gaia and our orbit analysis suggest that θ² Ori A is such a star. We also conclude that NU Ori was ejected from the ONC centre 0.1–0.2 Myr ago.

The dynamical ejection of stars can occur in any star cluster that contains multiple massive stars and is crucial for their formation because it changes the timing of the ionization of gas (Kroupa et al. 2018; Wang, Kroupa & Jerabkova 2019; Fujii et al. 2021c). In addition to the formation of a star cluster, ejected massive stars may affect star formation on a larger scale. The observed high runaway fractions of massive stars (20  per cent; Blaauw 1961) can be explained by the dynamical ejection from star clusters (Fujii & Portegies Zwart 2011). Runaway massive stars travel from their origin (dense molecular clouds) to low-density environments and efficiently ionize the interstellar gas. This may trigger new star formation (Andersson, Agertz & Renaud 2020). Some of these massive stars die in supernova explosions ∼100 pc (≃30 [km s⁻¹] × 3 [Myr]) away from their origin, which is a 1-pc-sized star cluster. This strong feedback can affect the evolution of the interstellar medium in the 100-pc scale, as has also been suggested for the Orion complex (Kounkel 2020), and even larger galactic scale star formation (Andersson et al. 2020).

ACKNOWLEDGEMENTS

The authors thank Steven Rieder for providing the amuse script for surface density and temperature map, Kurumi Ishikura and Ryoichi Nishi for discussion on the stellar distribution around the ONC, Takaaki Takeda (4D2U at the National Astronomical Observatory of Japan) for the visualization of the simulation, and Editage (www.editage.com) for English language editing. Numerical computations were carried out on Cray XC50 CPU-cluster at the Center for Computational Astrophysics (CfCA) of the National Astronomical Observatory of Japan. This work was supported by JSPS KAKENHI Grant Numbers 19H01933,20K14532,21J00153,21K03614,21K03633,21H04499 and Initiative on Promotion of Supercomputing for Young or Women Researchers, Information Technology Center, the University of Tokyo, and MEXT as ‘Program for Promoting Researches on the Supercomputer Fugaku’ (Toward a unified view of the universe: from large scale structures to planets, Revealing the formation history of the universe with large-scale simulations and astronomical big data). MF was supported by the University of Tokyo Excellent Young Researcher Program. LW thanks the support from the one-hundred-talent project of Sun Yat-sen University and the National Natural Science Foundation of China through grant 12073090. LW alsothanks the financial support from JSPS International Research Fellow (Graduate School of Science, the University of Tokyo).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author. petar is available here: https://github.com/lwang-astro/PeTar. CELib is available here: https://bitbucket.org/tsaitoh/celib.

APPENDIX A: OB STARS NEAR THE ONC FROM GAIA CATALOGUE

Table A1 lists the 34 OB stars within 5 pc (0|${_{.}^{\circ}}$|71) from the ONC.

Table A1.

Catalogue of OB stars in the ONC region.

Name	Sp. type	RA	Dec.	Parallax	(μα*, μδ)	(Δvα, Δvδ)	Δv
		(°)	(°)	(mas)	(mas yr−1)	(km s−1)	(km s−1)
Brun 508	B9V	+83.77	−5.98	2.62 ± 0.03	(+1.38, +0.53)	(+0.51, +0.42)	+0.65
HD 36919	B9V	+83.70	−6.00	2.64 ± 0.03	(+0.92, +0.07)	(−0.33, −0.42)	+0.53
iot Ori B	B8III	+83.86	−5.91	2.79 ± 0.05	(+1.13, +1.62)	(+0.05, +2.25)	+2.25
HD 37000	B3/5	+83.80	−5.93	2.62 ± 0.04	(+1.44, −0.07)	(+0.62, −0.67)	+0.92
HD 36983	B5(II/III)	+83.78	−5.87	2.63 ± 0.03	(−0.35, +0.59)	(−2.60, +0.53)	+2.66
HD 36999	B8(III)	+83.81	−5.83	2.60 ± 0.04	(+1.34, +0.52)	(+0.44, +0.39)	+0.59
HD 36917	B9III/IV	+83.70	−5.57	2.22 ± 0.06	(+2.77, −1.74)	(+3.57, −4.36)	+5.63
HD 36939	B7/8II	+83.73	−5.51	2.38 ± 0.04	(+1.04, +0.64)	(−0.12, +0.67)	+0.68
HD 37150	B3III/IV	+84.06	−5.65	2.66 ± 0.05	(+1.21, −0.15)	(+0.19, −0.80)	+0.82
HD 37174	B9V	+84.11	−5.41	2.63 ± 0.02	(+1.27, +0.46)	(+0.31, +0.29)	+0.42
V* V1073 Ori	B9.5V	+83.87	−5.44	2.61 ± 0.03	(−0.11, +1.01)	(−2.19, +1.29)	+2.54
HD 36982	B1.5Vp	+83.79	−5.46	2.45 ± 0.02	(+1.62, +1.78)	(+1.00, +2.87)	+3.04
tet02 Ori C	B4V	+83.88	−5.42	2.45 ± 0.04	(+2.51, +3.73)	(+2.73, +6.63)	+7.17
tet02 Ori B	B2-B5	+83.86	−5.42	2.39 ± 0.05	(+1.16, +0.16)	(+0.12, −0.27)	+0.30
tet02 Ori A	O9.5IVp	+83.85	−5.42	2.97 ± 0.21	(+1.09, +2.39)	(−0.02, +3.34)	+3.34
tet01 Ori D	B1.5Vp	+83.82	−5.39	2.28 ± 0.03	(+1.82, +0.39)	(+1.50, +0.19)	+1.51
tet01 Ori C	O7Vp	+83.82	−5.39	2.50 ± 0.14	(+2.26, +0.99)	(+2.20, +1.32)	+2.56
tet01 Ori A	B0V	+83.82	−5.39	2.64 ± 0.07	(+1.36, +0.25)	(+0.46, −0.09)	+0.47
V* V1230 Ori	B1	+83.84	−5.36	2.46 ± 0.03	(+3.06, −1.42)	(+3.79, −3.33)	+5.05
NU Ori	O9V	+83.88	−5.27	2.41 ± 0.06	(+0.92, +1.72)	(−0.35, +2.79)	+2.82
HD 36655	B9V	+83.28	−5.34	2.82 ± 0.04	(+0.11, +0.19)	(−1.67, −0.18)	+1.68
HD 36981	B7III/IV	+83.78	−5.20	2.59 ± 0.04	(+0.91, +0.37)	(−0.34, +0.13)	+0.36
HD 37060	(B9)	+83.89	−5.11	2.58 ± 0.02	(+1.40, +0.79)	(+0.54, +0.90)	+1.05
HD 37059	B8/A0V	+83.88	−4.90	2.60 ± 0.03	(+1.71, +0.81)	(+1.12, +0.93)	+1.45
HD 37058	B3/5II	+83.89	−4.84	2.62 ± 0.04	(+1.49, +0.74)	(+0.71, +0.80)	+1.07
HD 294264	B3	+83.81	−4.86	2.43 ± 0.05	(+1.42, −1.44)	(+0.63, −3.39)	+3.45
c Ori	B1V	+83.85	−4.84	6.52 ± 1.52	(−4.12, −1.97)	(−3.80, −1.65)	+4.14
HD 36938	B9V	+83.73	−4.77	2.51 ± 0.03	(+1.99, −0.71)	(+1.68, −1.90)	+2.54
HD 36958	B3/5V	+83.77	−4.73	2.76 ± 0.08	(+0.97, +0.11)	(−0.22, −0.32)	+0.39
HD 37130	B8/9IV	+84.01	−4.75	2.56 ± 0.02	(+0.22, −1.01)	(−1.63, −2.42)	+2.92
HD 37025	B3(III)	+83.82	−6.03	2.65 ± 0.08	(+1.39, +2.04)	(+0.52, +3.11)	+3.15
HD 36960	B1/2Ib/II	+83.76	−6.00	2.62 ± 0.12	(+1.11, +1.68)	(+0.02, +2.50)	+2.50
HD 36959	B1.2	+83.75	−6.01	2.79 ± 0.11	(+0.96, −1.00)	(−0.24, −2.21)	+2.23
HD 36918	B8.3	+83.70	−6.01	2.66 ± 0.03	(+1.20, +0.51)	(+0.19, +0.37)	+0.41

Note. From the left, columns 1 and 2 represent the star name and the spectral type in SIMBAD, respectively. Column 3 represents the right ascension. Column 4 represents the declination. Column 5 represents the parallax and its uncertainty in Gaia EDR3. Column 6 represents the proper motion vector in Gaia EDR3. Columns 7 and 8 represent the velocity vector with respect to the ONC’s centre of mass and its magnitude, respectively. HD 37025–HD 36918 are located within 5 pc from the ONC in projection but out of the region shown in panel (c) of Fig. 1.

In the following, we summarize the details of massive stars (>10 M_⊙) that are located in the Orion Nebula (M42 and M43). Here, we assumed the distance to the ONC and the radial velocity of the ONC as 388 ± 5 pc (Kounkel et al. 2017) and |$27.45^{+0.21}_{-0.22}$| km s⁻¹ (Theissen et al. 2022), respectively.

A1 θ¹ Ori C

θ¹ Ori C (HD 37022) is the most massive star in the ONC. The parallax from Gaia EDR3 is 2.50 ± 0.14, which is |$400^{+24}_{-21}$| pc (379–424 pc). The radial velocity is 24.50 ± 1.2 km s⁻¹ (Conti 1972). The proper motion relative to the ONC is 2.56 km s⁻¹. This value is similar to the velocity dispersion of stars in the ONC core measured with Gaia (Theissen et al. 2022). From these data, we consider that θ¹ Ori C is inside the ONC. In Fig. A1, we visualized the data of this star.

Figure A1.

Distribution of bright stars (G < 12 mag) near the ONC for which the stellar parallax data are consistent with 400 < d < 420 pc within the 3σ level. The blue circles and green squares represent stars classified in SIMBAD as O- and B-type stars, respectively. The magenta triangles represent other bright stars with relatively blue colour (G_BP − G_RP) < 0.5. The thick blue circle corresponds to θ¹ Ori C. Top left-hand panel: the relative proper motion vector of stars with respect to the ONC’s systemic proper motion in Jerabkova et al. (2019). Top right-hand panel: the colour–magnitude diagram. Bottom left-hand panel: the parallax and its uncertainty. Bottom right-hand panel: summary of the observational property of θ¹ Ori C and its tangential velocity in units of |$\rm {km\,\,s^{-1}}$| assuming a fiducial distance of |$d = 400 \,\,\rm {pc}$|⁠.

We calculated the mass-loss rate using equation (12) of Vink et al. (2000) with the stellar parameters summarized in Table 3 and obtained 4.0 × 10⁻⁷ M_⊙ s⁻¹. The obtained mechanical luminosity was 7.9 × 10³⁵ erg s⁻¹, which is consistent with previous research (7–8 × 10³⁵ erg s⁻¹; Güdel et al. 2008; Pabst et al. 2019).

A2 θ² Ori A

θ² Ori A (HD 37041) is the second brightest star in this region. The estimated mass is 25–39 M_⊙ (Preibisch et al. 1999; Simón-Díaz et al. 2006). The parallax is 2.97 ± 0.21, which corresponds to the distance of |$336\pm ^{+26}_{-22}$| pc (314–362 pc). This star is located in front of the ONC, albeit the distance is relatively uncertain (see Fig. A2). The radial velocity is 20–30 km s⁻¹ (Henroteau & Henderson 1920; Aikman & Goldberg 1974).

Figure A2.

Same as Fig. A1 but for θ² Ori A.

From the distance including the error, this star is at least 20 pc closer to observer with 1σ error. However, the parallax error of this star is ∼10  per cent, and the error of 0.1 mas in parallax corresponds to 16 pc at a distance of 400 pc. Therefore, the distance of θ² Ori A to the ONC can be shorter. The projected distance of θ² Ori A to the ONC is ∼0.3 pc. The projected relative velocity is 3.34 km s⁻¹, and the direction of the velocity is towards the ONC centre.

We also estimated the mechanical luminosity of θ² Ori A. The mass-loss rate is obtained to be 5.9 × 10⁻⁸ M_⊙ yr⁻¹ with the stellar parameters summarized in Table 3. We assume the wind terminal velocity as 2000 km s⁻¹ (Prinja et al. 1990). The estimated mechanical luminosity is 7.4 × 10³⁴ erg s⁻¹, which is an order of magnitude smaller than that of θ¹ Ori C.

A3 NU Ori

NU Ori (HD 37061) is an O9V- or B0.5V-type star at the centre of the small H ii region in M43 (Simón-Díaz et al. 2011). The parallax is 2.41 ± 0.06. The corresponding distance is |$415^{+11}_{-10}$| pc (405–426 pc). The observed radial velocity is 66.90 km s⁻¹ (SIMBAD; Abt 1970), which suggests >30 km s⁻¹ to the ONC. The proper motion relative to the ONC is 2.82 km s⁻¹, and it is apparently escaping from the ONC (see Figs 1 and A3). This suggests that this star may be a runaway star returning towards the ONC.

Figure A3.

Same as Fig. A1 but for NU Ori (HD 37061).

The mass is estimated to be 13–16 M_⊙ (Gravity Collaboration et al. 2018). From the stellar parameters given in Table 3, we obtained L_w = 3.1 × 10³³ erg s⁻¹.

APPENDIX B: SIMULATION WITHOUT GRAVITATIONAL SOFTENING

We also confirmed that the dynamical ejections of massive stars are crucial for the formation of star clusters by performing a reference run with a gravitational softening, with which the strong dynamical encounters among stars are ‘softened’ and as a result the ejection of stars is suppressed.

We additionally performed a simulation with gravitational softening for stars. In this simulation, the initial condition of gas is exactly the same as Model A. The softening length for stars was set to be the same as that for gas particles (0.07 pc). The snapshots at t_sim = 6.5 Myr are shown in Fig. B1. Unlike in the case without softening (see Fig. 2), some small stellar clumps begin ionizing the surrounding gas from inside the clumps. This prevents them from evolving into a massive cluster via mergers of the clumps. With a gravitational softening, massive stars born in small clumps are kept inside the clumps without being ejected to the outside of the clumps. Therefore, the feedback inside the clumps works more efficiently when compared with the case without softening. We observed the suppression of stellar feedback owing to the ejection of massive stars. This was discussed in Kroupa et al. (2018).

Figure B1.

Snapshots of the simulation with a softening. The softening length is 0.07 pc, which is the same as that for gas particles. Time indicates the time from the beginning of the simulation.

Author notes

JSPS Research Fellow.