Turbulence in Milky Way Star-forming Regions Traced by Young Stars and Gas

Trung Ha; Yuan Li; Marina Kounkel; Siyao Xu; Hui Li; Yong Zheng

doi:10.3847/1538-4357/ac76bf

1. Introduction

The interstellar medium (ISM) is multiphase and turbulent. Interstellar turbulence affects the velocity statistics on all scales from the diffuse ionized gas to the dense molecular gas. Larson (1981) compiled the internal velocity dispersion of nearby molecular clouds from the width of various emission lines and found a power-law relation between the velocity dispersion and cloud size, with a slope of 0.38, between the 1/3 slope expected for turbulence dominated by solenoidal motions (Kolmogorov1941) and the 1/2 slope expected for shock-dominated turbulence (e.g., Galtier & Banerjee 2011; Federrath 2013).

Inside molecular clouds, massive O- and B-type stars create H ii regions by ionizing and heating the ISM with UV radiation (see, e.g., Spitzer 1978). Compression from supernova (SN) explosions can trigger further star formation (Kennicutt & Evans 2012). Interstellar processes such as SNe, stellar winds, ionizing radiation, and cosmic rays all play a role in maintaining the structure of the ISM, as well as in regulating star formation (Hennebelle & Falgarone 2012; Rathjen et al. 2021; Xu & Lazarian 2022; Hu et al. 2022). Of these mechanisms, SN explosions are perhaps the most significant source of energy injection into the local ISM (Mac Low & Klessen 2004; Padoan et al. 2016).

Many numerical simulations have been carried out to study how SN explosions affect the surrounding ISM and star formation in the last decade. In general, they find that most of the kinetic energy from an SN explosion is transferred to the low-density region of the ISM in both single- (see, e.g., Kim & Ostriker 2015; Li et al. 2015; Martizzi et al. 2015) and multiple-SN events (Padoan et al. 2016), suggesting that turbulence may not be evenly driven in different phases of the ISM.

Observational studies of ISM turbulence have been done with various gas tracers. Measurements of turbulence through the density spectrum were carried out using the Wisconsin H-Alpha Mapper (WHAM; Haffner et al. 2003, 2010), which revealed a turbulent cascade consistent with Kolmogorov theory (Chepurnov & Lazarian 2010). Burkhart et al. (2010) found that most of the H i in the Small Magellanic Cloud (SMC) is subsonic or transonic and concluded that turbulence in the SMC is dominated by hydrodynamical pressure over magnetic pressure.

Statistical measurements of velocity fluctuations in the neutral ISM can be done using the Velocity Channel Analysis, the Velocity Coordinate Spectrum (Lazarian & Pogosyan 2000, 2006), the Delta-variance technique (Stutzki et al. 1998; Ossenkopf & Mac Low 2002), or the Principal Component Analysis (Heyer & Brunt 2004; Roman-Duval et al. 2011) based on spectroscopic data. These analyses generally reveal power-law spectra of turbulent velocities (Lazarian 2009; Chepurnov et al. 2010).

Additionally, the kinematics of Galactic H ii regions have been studied using the observations of ionized gas. Fabry–Perot interferometry was used to probe the Hα line in emission nebulae. Louise & Monnet (1970) found a velocity–size relation that agrees with the 1/3 slope of Kolmogorov turbulence in M8, while Roy & Joncas (1985) and Joncas & Roy (1986) found slopes of the structure functions closer to 1/2, indicating turbulence in the supersonic regime in both S142 and M17, respectively.

Closer to our solar neighborhood, studies into the nature of turbulence in the Orion Nebula have been an ongoing effort since the 1950s. von Hoerner (1951) observed an rms difference in emission-line velocities that increases with larger angular separations. He reasoned that such a relation could be explained by the Kolmogorov theory if the effective depth of the nebula was sufficiently small as to ignore projection effects. Subsequent studies, carried out by Münch (1958) and Castaneda (1988), dispute this conclusion. Münch (1958) found an optical depth smaller than the value postulated by von Hoerner (1951) and attributed the failure of Kolmogorov theory to the compressibility of turbulence in the Orion Nebula, while Castaneda (1988) argued against a simple power-law relation, instead suggesting energy injections at different scales.

ISM turbulence has also been probed using "point sources" such as molecular cores. Qian et al. (2012) developed statistical measurements of core-to-core velocity dispersion over a range of length scales, i.e., the core velocity dispersion technique. Qian et al. (2018) find signatures of a turbulent power-law spectrum in the Taurus Molecular Cloud and demonstrate that the statistics of turbulent velocities are imprinted in the velocities of dense cores (Xu 2020).

One limitation of using gas tracers is they lack three-dimensional information on the position, and the velocity information is also limited to only the line of sight. In Ha et al. (2021, hereafter Paper I), we proposed a new method to measure turbulence in the ISM with young stars. Because stars are formed out of turbulent molecular clouds, we expected these young stars to inherit the turbulent kinematics from their natal clouds. Using six-dimensional (6D) position and velocity information of >1400 stars in the Orion Complex, we found that the first-order velocity structure function (VSF) of young stars in loose groups exhibits characteristics of turbulence. Our method can therefore be used to probe the turbulent properties of molecular clouds in addition to the analysis of gas kinematics and has been used in recent work by Krolikowski et al. (2021), Zhou et al. (2022).

In this work, we apply a similar analysis to a larger sample of molecular complexes, which include three of the nearest massive star-forming regions: Ophiuchus, Perseus, and Taurus. Additionally, we analyze the kinematics of Hα gas from the WHAM survey to study the connection between stars and the multiphase ISM.

In Section 2, we describe the data used in this analysis, the reduction pipeline for the Hα kinematic survey, and how we compute the VSFs using stars and Hα. In Section 3, we present the VSF of young stars and the multiphase gas of all four regions. In Section 4, we discuss sources of uncertainty associated with our analysis and the nature of turbulence in the Milky Way ISM. We summarize our results in Section 5.

2. Observations, Data Reduction, and Methods

2.1. Observations and Data Analysis

2.1.1. Stellar Observations

We obtain the locations, parallax distances, and proper motions of stars from the Gaia mission (Gaia Collaboration et al. 2016) in four nearby star-forming regions: the Orion Molecular Cloud complex (Orion), the ρ Oph dark cloud (Ophiuchus), the Perseus molecular cloud traced by two clusters, NGC 1333 and IC 348 (Perseus), and the Taurus molecular clouds (Taurus). The Gaia catalog includes five-parameter astrometric measurements (plane-of-sky coordinates, proper motions, and parallax distances) of over 1.46 billion stars in the Milky Way. Only a very limited subset (∼7 million) has line-of-sight velocities measured by Gaia. The five-dimensional astrometric data for our four regions are obtained from Gaia's third early data release (EDR3; Gaia Collaboration et al. 2021).

The line-of-sight velocities of the stars in our study were observed in the near-infrared with the Apache Point Observatory Galactic Evolution Experiment (APOGEE) spectrograph, mounted on the 2.5 m Sloan Foundation Telescope of the Sloan Digital Sky Survey. APOGEE-2 Data Release 17 (DR17) collected high-resolution spectral images of over 657,000 stars in the Milky Way (Gunn et al. 2006; Blanton et al. 2017; Abdurro'uf & Aerts 2022). Combining these two surveys provides 6D information (3D positions and 3D velocities) of stars for the four molecular clouds.

The stars in our study were assigned to each of the four regions by Kounkel & Covey (2019). They performed a clustering analysis on Gaia DR2 sources in the Milky Way using Python implementation of HDBSCAN (Hierarchical Density-Based Spatial Clustering of Applications with Noise; McInnes et al. 2017). We crossmatch this catalog with Gaia EDR3 and APOGEE-2 to obtain 6D astrometries of the stars.

Then, we exclude a small number of stars whose parallax distances d, line-of-sight velocities v_LSR, and line-of-sight velocity errors v_LSR,err deviate by more than 4σ from the mean of the distribution. In Taurus, one star with very high reported proper motion (v_LSR,RA > 100 mas yr⁻¹) is excluded. The ranges of d and v_LSR, and the threshold of v_{LSR, err} that we use, as well as the total number of stars before and after the cut, are specified in Table 1. We note, however, that the exact cut of these ranges does not affect the conclusion of this work.

2.1.2. Hα Observations

The Hα data in our study come from the WHAM survey (Haffner et al. 2003, 2010). We obtained the data cubes from the "DR1-v161116-170912" release, which comprises more than 49,000 spectra divided into an irregular grid with a spatial resolution of 1°, and the data cube is Nyquist-sampled with 4 × 4 pixels per resolution element (Haffner et al. 2003). The intensity of the Hα emission is recorded in a three-dimensional array encoding its Galactic longitude, Galactic latitude, and line-of-sight velocity (−170 km s⁻¹ < v_LSR < 160 km s⁻¹). We fit a Gaussian distribution to the flux profile along individual lines of sight. Then, we take the velocity at the peak of the distribution as the centroid velocity and the uncertainty of the fit as the uncertainty. We also correct for Galactic rotation following the procedure in Qu et al. (2020).

Figure 1 shows the Hα map of the four regions in our analysis, with stars overplotted as foreground dots. The top panels show the integrated line-of-sight fluxes in log scale. In the foreground, each star is color-coded by its distance from Earth. We note that although Perseus and Taurus are very close to each other on the sky, they are physically separated by ∼150 pc along the line of sight. In the bottom panels, we plot the centroid velocity of the Hα in the background, and the line-of-sight velocities of stars in the foreground.

2.2. Data Analysis

We compute the first-order VSF of stars in each region. The VSF is related to the kinetic energy power spectrum of a velocity field and is expressed as a two-point correlation function between the absolute velocity differences ∣δ v∣ versus the physical separation ℓ. For a turbulent velocity field, we expect the VSF to have a characteristic power-law scaling with a slope of 1/3 for incompressible turbulence dominated by solenoidal motions (Kolmogorov 1941) or a slope of 1/2 for supersonic turbulence dominated by compressive motions. For the Milky Way molecular clouds, Larson (1981) found a relationship between cloud size and velocity dispersion to be σ = 1.1 R^0.38.

The 3D VSF of stars is computed as follows. First, the positions, distances, radial velocities, and proper motions of stars are translated into galactocentric Cartesian coordinates and velocities. Then, for each pair of stars, the absolute velocity vector difference ∣δ v ∣ and the distance d between them is calculated. For individual clouds, these distances are sorted into logarithmic bins of ℓ, and ∣δ v ∣ values within a bin are averaged out to obtain 〈∣δ v ∣〉. To estimate uncertainties, we first generate 1000 random samples of stellar parameters following a Gaussian distribution around the reported measurements and with a 1σ scatter equal to the reported uncertainties. Then, at each ℓ, we take the 1σ value of the 1000 〈∣δ v ∣〉 to be the uncertainty.

For the WHAM Hα data, we first define a rectangular box containing the stars in each of the four regions as specified in Table 1 and shown in Figure 1. Our results are not sensitive to the exact size of the box (see Appendix C for detailed discussions). We compute the first-order VSF of each region using the projected position–position–velocity information. The resultant VSF is 〈∣δ v∣〉_proj versus angular separations θ. Then, we assume a constant distance equal to the mean distance of stars in each region (see Table 1) and translate θ (°) to ℓ_proj (pc). Because three degrees of freedom are missing in this projected VSF calculation (the proper motions of gas and its exact distance), we scale 〈∣δ v∣〉_proj and ℓ_proj by a constant factor, assuming an isotropic turbulent velocity field:

$\begin{eqnarray}&&\left\langle \left|\delta {\boldsymbol{v}}\right|\right\rangle =\sqrt{3}\,{\left\langle \left|\delta v\right|\right\rangle }_{\mathrm{proj}},\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&{\ell }=\sqrt{\displaystyle \frac{3}{2}}\,{{\ell }}_{\mathrm{proj}}.\end{eqnarray} \tag{ 2 }$

**Figure 1.** Top panels: positions of young stars color-coded by their distances in Orion, Ophiuchus, and Taurus-Perseus (from left to right). The background shows the Hα flux from the WHAM survey (lighter is brighter). The orange and green boxes denote the areas where Hα is analyzed for Taurus and Perseus, respectively. The bottom panels show the line-of-sight velocities of the young stars (in the LSR frame) in these regions. The background shows the Hα line-of-sight centroid velocities.
Download figure:
Standard image High-resolution image

Table 1. Properties of Star-forming Regions

Region	Number	Number of	Average	GLON	GLAT	Distance	v_LSR Range	v_LSR,err Cut
	of Stars	Stars Used	Distance (pc)	(°)	(°)	Range (pc)	(km s⁻¹)	(km s⁻¹)
Orion	2647	2468	410	190, 217	−28, −5	300, 500	−10, 40	10
Ophiuchus	113	107	140	342, 360	12, 28	120, 170	−30, 30	2
Perseus	159	138	320	156, 164	−23, −15	250, 400	−6, 30	10
Taurus	157	139	140	165, 177	−23, −11	97, 200	−5, 25	5

Note. The number of stars is all the stars with 6D information. The number of stars used is the stars in our analysis after excluding outliers. GLON (Galactic Longitude) and GLAT (Galactic Latitude) ranges denote the borders of windows used to compute the Hα VSF. The distance ranges, v_LSR ranges, and v_{LSR, err} cuts are applied to the stars to exclude outliers.

Download table as: ASCII Typeset image

3. Results

3.1. The Universality of ISM Turbulence Traced by Young Stars

We first examine the velocity statistics of the four star-forming regions: Orion, Ophiuchus, Perseus, and Taurus. Figure 2 shows the VSFs traced by stars of these four regions. For reference, we also plot the Larson's relation of velocity dispersion versus cloud size, the 1/3 slope of the Kolmogorov turbulence, and the 1/2 slope of the supersonic turbulence. Immediately, we recognize a remarkable universality among the VSFs of the different regions. Over a large dynamical range, the magnitude and slopes of the VSFs show a general agreement with Larson's relation (see Table 2 for the best-fit slopes of the VSFs). This suggests that the motion of young stars in the Milky Way shows characteristics of turbulence and supports our previous findings in Paper I that turbulence in the ISM can be traced with young stars.

Table 2. Slopes of the VSFs

Region	Slope of	Fitting	Slope of the	Fitting
	Stars' VSF	Range (pc)	Hα VSF	Range (pc)
Orion	0.25 ± 0.01	12−75	0.55 ± 0.02	9−150
Ophiuchus	0.34 ± 0.03	7−28	0.36 ± 0.03	3−30
Perseus	0.22 ± 0.03	5−80	0.61 ± 0.08	7−55
Taurus	0.34 ± 0.03	7−70	0.42 ± 0.12	3−25

Note. Inside the fitting ranges, the VSFs are fitted with an equation: 〈∣δ v ∣〉 = A · ℓ^B, where B is the slope. The uncertainties reported with each slope are obtained from the covariance matrix, taking into account the uncertainties in 〈∣δ v ∣〉 shown in Figure 3. We choose the fitting ranges such that the uncertainties are low and there is a consistent trend in the slope of the VSF.

Download table as: ASCII Typeset image

At ℓ ≲ 10 pc, the slopes of the VSF of all four regions are flatter than those of Larson's relation, as well as the Kolmogorov (slope 1/3) turbulence. We discuss the possible cause of this flattening in Section 4.1.

3.2. The VSFs of Individual Regions

Figure 3 shows the VSFs of young stars as well as gas in each of the four regions: Orion, Ophiuchus, Perseus, and Taurus. The best-fit slopes of these VSFs are listed in Table 2.

3.2.1. Orion

The first-order VSF of Orion traced by stars has been studied in detail in Paper I. We only reiterate some of the main points here and address the differences caused by the addition of ∼1000 stars compared with Paper I. Throughout the inertial range, the magnitude of $\langle | \delta \vec{v}| \rangle$ is higher than Larson's relation. At small scales (1 pc ≲ ℓ≲ 20 pc), the number statistics is dominated by stars in the Orion Nebula Cluster (ONC), which flattens the VSF (see Section 4.1 for more discussion on the small-scale flattening). At ranges from 20 to ∼70 pc, Orion D, which is a loose stellar association at the front of the Orion Complex along our line of sight, dominates the statistics, and it exhibits a power-law scaling, which indicates an imprint of turbulence from their natal molecular clouds. The VSF of Orion without considering the ONC would be steeper at small scales, which renders it more in line with the VSF of other regions (see Figure 2). At ∼70 pc, there is a bump in the structure function, and we proposed that this is evidence supporting the existence of an SN explosion at the heart of Orion 6 Myr ago (Kounkel 2020). Paper I used Gaia DR2 only. The newly added stars in this study from Gaia EDR3 have larger overall uncertainty in their parametric solutions. Therefore, the noise floor for Orion at small scales is boosted from ∼3 km s⁻¹ to ∼5 km s⁻¹. This slightly elevates the overall VSF compared to Paper I.

Hα emission in Orion mainly comes from Barnard's loop and the λ Ori bubble (see Figure 1), which are thought to be created by SN explosions (Kounkel 2020). We see the effect of such energy injection in the VSF traced by Hα. At ℓ ≳ 20 pc, the amplitude of the Hα VSF is higher than that of the VSF of young stars. Over the whole dynamical range we probe here, the Hα VSF has a roughly constant slope of 1/2.

The VSF traced by Hα shows a bump at ℓ ∼ 150 pc. This is indicative of local energy injection, similar to those seen in the VSF traced by young stars. However, we caution that this scale is very close to the size of the Hα box, and the bump can also be an artifact (see, e.g., Mohapatra et al. 2022). We discuss this effect in more detail in Appendix C. At very small scales (θ = 1°, corresponding to ℓ ∼ 7 pc), our analysis is limited by the resolution of the WHAM survey (Haffner et al. 2003). We discuss the effect of the resolution limit and other uncertainties related to the Hα VSF in Section 4.2.

3.2.2. Ophiuchus

The young stars in Ophiuchus have a VSF similar to that of the Orion Complex. At small ℓ, the structure function is flatter than 1/3 and then steepens to ∼1/3 beyond 10 pc up to ∼25 pc. The Hα VSF in Ophiuchus has a slope close to 1/3 and a higher amplitude than the VSF of the stars. Both VSFs show clear bumps at ∼25 pc, which we interpret as evidence of local energy injection.

In Orion, such an energy injection is associated with a past SN explosion as we discuss in detail in Paper I. If the energy injection in Ophiuchus is also caused by an SN, we estimate that the event took place around 25 pc/v_region ∼ 2.5 Myr ago, where we use the characteristic Hα velocity v_region ∼ 10 km s⁻¹. Such an estimate is broadly consistent with observations of runaway stars that trace back to an SN explosion in the region ∼2 Myr ago (Neuhäuser et al. 2020).

We note that the scale of the bump is close to the size of the region containing the young stars in our analysis. Therefore, it is technically possible that the bump in the VSF of young stars is only an artifact caused by the limited sample size. In Appendix A, we analyze the entire Upper Sco region. We find that as long as we apply our analysis to young stars (<2 Myr), the results are consistent, while the VSFs of older stars do not show this feature. In addition, the SN scenario is also supported by the bump in the Hα VSF.

The cold and dense molecular cores, where stars are born out of, emit CO. Qian et al. (2015, 2018) applied the core velocity dispersion (CVD; Qian et al. 2012) method to analyze turbulence in star-forming regions. We translate their best-fit CVD to a VSF of CO by $\mathrm{VSF}=\sqrt{\tfrac{2}{\pi }}\,\mathrm{CVD}$ .

In Ophiuchus, the amplitude of the CO VSF (pink solid lines in Figure 3) is much lower than the VSF of young stars and Hα. It is possible that the recent energy injection has affected the warm phase of the ISM more than the cold phase (see discussion in Section 4.3). It is also possible that the low amplitude is due to the spatial variation of ISM turbulence (see discussion in Appendix C). Qian et al. (2015) only examined the smaller-scale velocity structures (ℓ≲ 10 pc), where there is no reliable component that we can compare against (due to the relaxation of stars and the resolution limit of Hα).

3.2.3. Perseus

At length scales below 10 pc, the VSF of young stars in Perseus is very similar to that of Ophiuchus (see Figure 2). Because we include two clusters in our analysis of Perseus: NGC 1333 and IC 348, the VSF is flattened from internal dynamical relaxation. The existence of subgroups and subclusters is likely the cause of the small wiggles in the VSF. At larger scales, the VSF is generally consistent with Larson's relation.

The VSF traced by Hα here has an amplitude similar to the VSF traced by young stars but with a steeper slope (∼1/2). VSF traced by CO matches the Hα VSF at ℓ > 6 pc but flattens on smaller scales due to projection effects for thick clouds (Qian et al. 2015; Xu 2020).

3.2.4. Taurus

Taurus is overall very similar to Perseus. The VSFs of young stars, Hα, and CO generally agree with each other, as well as Larson's relation.

There is no prominent bump detected in the VSF in both Perseus and Taurus, which is consistent with the findings of Krolikowski et al. (2021) that there have been no SN activities in the last 5 Myr in these regions. However, there are wiggles scattered throughout the length scales, likely because Taurus contains several small dense groups that are primarily homogeneously distributed throughout the region.

4. Discussion

4.1. Uncertainties and Biases in Turbulence Traced by Young Stars

We discussed various biases and uncertainties associated with measuring turbulence from stars in Paper I, including dynamical relaxation, binary contamination, and drifting. In this section, we reiterate the main points and discuss additional mechanisms that can affect the shape of the VSF.

Stellar interactions can cause dynamical relaxation, which erases the memory of turbulence and flattens the VSF. Similar to Orion's loose groups discussed in Paper I, the relaxation times of Ophiuchus, Perseus, and Taurus are estimated to be tens of Myr, much longer than their age (for a detailed discussion of the age of each region, see Wilking et al. 2008; Pavlidou et al. 2021; Krolikowski et al. 2021). Thus, dynamical relaxation has not affected the overall shape of the VSF for the three regions. However, there are some small, compact star clusters in these regions, which can contribute to the flattening at small scales (ℓ ≲ 10 pc; see Figure 2).

Another source of uncertainty originates from the drifting of stars along their natural trajectories over time. Stars born with high velocity differences drift apart at faster rates than stars born with low velocity differences, which steepens the VSF roughly over the crossing time (t_cross) of the cloud. Similar to the star groups in Orion that we discussed in Paper I, the t_cross of Ophiuchus, Perseus, and Taurus are all larger than their current age, and thus, this effect is unlikely to have significantly altered the shape of the structure functions.

Low-velocity binary systems and contamination from field stars can act like noise in our analysis, which also tends to flatten the VSF as discussed in Paper I.

Additionally, other stellar feedback processes such as winds, jets, and radiation can also inject energy on small scales (e.g., Gritschneder et al. 2009; Offner & Arce 2015; Li et al. 2019; Rosen et al. 2021; Grudić et al. 2022). This can contribute to the flattening of the VSFs at small scales, provided that these stars are massive (>8 M_⊙) and within ∼1 pc of each other (Rosen et al. 2021).

Stutz & Gould (2016) quantified an ejection mechanism through which protostars decouple from the molecular gas. They found that the undulation of star-forming filaments in Orion A "slingshots" protostars out of the filament at ∼2.5 km s⁻¹, which is also at the level where the VSF flattens in Ophiuchus, Perseus, and Taurus.

All sources of uncertainties listed above can contribute to the flattening of the VSF at small scales (ℓ ≲ 10 pc). Therefore, we only consider young stars to be reliable tracers of ISM turbulence above ∼10 pc.

The completeness of the stellar sample may also affect the VSFs. In Orion, the current stellar population with 6D measurements is estimated to cover about 25% of the total number of young stars in the region (Kounkel & Covey 2019). In the other regions, the coverage of our sample is ∼20%–30% (see Luhman et al. 2016; Krolikowski et al. 2021; Briceño-Morales & 2022).

If more young stars are observed in the future with levels of uncertainties similar to the existing data, we expect better statistics. Therefore, having a more complete sample of young stars in these regions can further strengthen the conclusion of our analysis. However, having a more complete sample of the young stellar population does not necessarily translate to better accuracy of the VSF analysis. In Orion, we obtained an additional ∼1000 stars going from Gaia DR2 (Paper I) to Gaia EDR3 (this work). However, these new stars have higher astrometric uncertainties compared to the ones previously identified in Gaia DR2, and thus the noise level of the VSF is elevated, and the slope of the VSF flattens compared to what we found in Paper I.

4.2. ISM Turbulence Traced by Hα

In Figures 2 and 3, we find that while the VSFs of stars are broadly consistent across the regions, the VSFs of Hα are more diverse.

In regions with recent SN explosions (Orion and Ophiuchus), the VSF of Hα gas has a higher amplitude than that of stars. Whereas in Perseus and Taurus, the Hα flux is lower (Figure 1) and the Hα VSF is more in line with the VSFs of stars and CO.

Within each individual region, we divide the pixels in half based on the Hα brightness—the high-emission half with all the pixels brighter than the median and the low-emission half with all the pixels dimmer than the median. Figure 4 shows the VSFs of high-emission Hα (dotted–dashed line) and low-emission Hα (solid line) within each region. In general, the VSF of high-emission Hα has a higher amplitude compared to that of the low-emission Hα.

**Figure 4.** VSF of four regions traced by Hα. Solid lines denote the VSF computed for pixels with low emission, and dotted–dashed lines denote the VSF computed for pixels with high emission.
Download figure:
Standard image High-resolution image

One possible explanation for this is that bright Hα comes from photoionization by young massive stars (Tacchella et al. 2022). Therefore, it traces the most intense stellar feedback regions and shows a higher VSF amplitude, whereas the low-emission Hα comes primarily from collisional ionization. Note that Orion has its high- and low-emission Hα VSFs essentially overlapping for the most part. This is likely because star formation in Orion is the most intense among all the regions studied here (see Table 2 of Lada et al. 2010), and thus most of the Hα is photoionized by massive stars.

There are many uncertainties and potential biases associated with the Hα analysis. One source of uncertainty comes from the resolution limit. Although the WHAM data cubes are sampled every 0 fdg 25, the reported spatial resolution of the survey is 1° (Haffner et al. 2003). At separations θ ≤ 1° (or ℓ_proj ≤ d · 1°), we treat the VSF traced by Hα as unresolved (represented by the unconnected dots in Figure 3). We note that the resolution limit might not significantly alter the overall VSF, as suggested by Li et al. (2020). Future studies with mock simulated Hα observations may help properly evaluate the effects of the resolution limit of the data set.

There are also uncertainties with the distance to the Hα gas. The bright Hα emissions in Orion and Ophiuchus are likely associated with the star-forming regions. However, because Perseus and Taurus do not emit intense Hα (see top right panel of Figure 1), it is uncertain whether the centroid velocities found in these two regions are at the exact same distance as the stars.

Another bias can come from the projection effects. When analyzing the VSF of Hα, we assume that all the gas is located at the same distance, effectively treating it as a 2D structure. This can flatten the VSF for thick clouds, as discussed in detail in Qian et al. (2015), Xu (2020), and Mohapatra et al. (2022). All the Hα VSFs in our analysis are either as steep as 1/3 or steeper. This suggests that either most of the Hα comes from a thin sheet of gas or the true VSF has an even steeper slope.

Another potential bias related to projection is that we are trying to infer 3D turbulence using projected 2D information. This can introduce uncertainties as the ISM is magnetized. In the presence of magnetic fields, turbulence is no longer isotropic (Seifried & Walch 2015; Hu et al. 2021).

We find that our results are not sensitive to the rotation model (see detailed discussions in Appendix B).

4.3. The Complex Nature of the ISM Turbulence

The ISM is multiphase in nature. The kinematic coupling between these phases is not well understood (for a review of the multiphase ISM model, see Cox 2005). More specifically, it is suggested that each phase possesses varying levels of turbulence, where ionized gas has been observed to have higher velocity dispersion than molecular gas in both high- and low-redshift galaxies (Ejdetjärn et al. 2022; Girard et al. 2021). In our study, we observe a similar trend comparing the VSF of Hα and other tracers of turbulence.

In Orion and Ophiuchus, there is a large deviation between the level of turbulence seen by different tracers. The VSFs traced by Hα show the largest amplitude at large scales, followed by the VSFs traced by stars, and finally by CO. One possible explanation is the difference in the transport of energy and momentum from SNe into different phases in these regions, as mentioned in Section 3.2. In an inhomogeneous medium such as the Milky Way ISM, SNe carry shocked matter radially outward following the path of least resistance. As a consequence, the cool dense phase traced by CO receives less energy injection than the warm diffuse phase traced by Hα, resulting in lower levels of turbulence in the corresponding VSF. Such behavior is also seen in hydrodynamical simulations of single SN remnants (e.g., Li et al. 2015; Martizzi et al. 2015).

In Perseus and Taurus, the VSFs traced by Hα and CO show a good agreement in both the amplitudes and slopes in these regions, except for the range ℓ ≲ 7 pc in Perseus, which suffers from projection effect due to the thickness of the cloud (see Section 4.2 of Qian et al. 2015, for a detailed discussion). At larger ℓ, the VSF traced by Hα and stars also closely follow each other. Such a trend is consistent with the evolutionary history of these regions with no clear imprint of recent SN explosions. Without this important local driver of turbulence, different phases of gas are well coupled, and the young stars formed out of these molecular clouds also largely retain the same memory of turbulence.

At length scales beyond the bumps caused by SN energy injection, we recognize a common upward trend in the VSFs of stars. In all four regions, the VSFs at the largest ℓ maintain a positive slope with no apparent driving scale. At these length scales, the ISM turbulence is driven primarily by galactic shear (see, e.g., Gammie et al. 1991) and accretion flows onto the galaxy (Forbes et al. 2022).

The difference in Hα VSFs in different regions shows that ISM turbulence can be unevenly distributed spatially as a result of local drivers. Turbulence can also be unevenly distributed between different phases of the ISM (e.g., Orion and Ophiuchus). Section 4.2 and Figure 4 further demonstrate that even with the same gas tracer Hα, the level of turbulence can vary depending on its brightness. In addition, we note that local energy injection from SN explosions is also intermittent temporally. Using the Hα VSF in Ophiuchus, we estimate the eddy turnover time at the driving scale to be

$\begin{eqnarray}&&\frac{\ell }{v}\sim \frac{20-30\,\mathrm{pc}}{5\,\mathrm{km}\,{{\rm{s}}}^{-}1}\sim 4-6\,\mathrm{Myr}.\end{eqnarray} \tag{ 3 }$

The Hα VSF in Orion gives a similar estimate. Typically, it takes a few eddy turnover times to establish a steady-state turbulent flow. On the other hand, given the size of these regions, the average time interval between SNe is on the order of a few Myr (Narayan & Ostriker 1990; Li et al. 2015). Each event also injects energy as a pulse with a very short duration. Therefore, both the "on" and "off" states of the driver are short compared with the timescale required to achieve a steady state. Thus, ISM turbulence in molecular complexes is often not in a steady state.

5. Conclusions

In this work, we analyze the kinematics of young stars and gas in four nearby star-forming regions: Orion, Ophiuchus, Perseus, and Taurus. For each region, we compute the first-order VSF of young stars and Hα gas and compare the results against each other, as well as analysis of CO molecular cores from the literature.

We find that the VSFs traced by young stars in all four regions exhibit similar properties, largely in agreement with Larson's relation, and reveal a universal scaling of turbulence in the ISM. This also confirms our previous finding in Ha et al. (2021) that young stars retain the memory of turbulence in their natal molecular clouds and thus can be used to measure the turbulent kinematics of molecular clouds.

The VSFs of Hα are more diverse, likely as a result of local energy injection. In regions with recent SN activities such as Orion and Ophiuchus, the VSFs of Hα are higher in amplitude compared to those without. The Hα VSFs are also higher than the VSFs of stars and CO in these regions, while regions without recent SN explosions (Perseus and Taurus) show well-coupled turbulence between different phases. Within each individual region, the high-emission Hα tends to show higher VSFs than the low-emission Hα.

Our findings support a complex picture of the Milky Way ISM, where turbulence can be driven at different scales and inject energy unevenly into different phases. Future observations with more tracers of turbulence and on more star-forming regions will help us better understand the diversity and complexity of ISM turbulence. In addition, high-resolution numerical studies can help further our understanding of the coupling between different phases and relevant timescales.

We thank Shmuel Bialy, John Forbes, Di Li, Mordecai-Mark Mac Low, Lei Qian, Amiel Sternberg, and Yuan-Sen Ting for the helpful discussions. We thank Zhijie Qu for sharing the Galactic rotation model. We also thank the anonymous referee for the helpful suggestions. H.L. and S.X. are supported by NASA through NASA Hubble Fellowship grants HST-HF2-51438.001-A and HST-HF2-51473.001-A, respectively, awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. Y.L. acknowledges financial support from NSF grant AST-2107735 and the College of Science and College of Engineering at UNT. Y.Z. is supported by NASA through a grant (HST-AR-16640.001-A) from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. The Wisconsin Hα Mapper and its Hα Sky Survey have been funded primarily by the National Science Foundation. The facility was designed and built with the help of the University of Wisconsin Graduate School, Physical Sciences Lab, and Space Astronomy Lab. NOAO staff at Kitt Peak and Cerro Tololo provided on-site support for its remote operation.

Facilities: Gaia - , Sloan - , WHAM. -

Software: astropy (Astropy Collaboration et al. 2013, 2018), matplotlib (Hunter 2007), numpy (van der Walt et al. 2011; Harris et al. 2020), pandas (McKinney 2010), scipy (Virtanen et al. 2020), Sagitta (McBride et al. 2021).

Appendix A: Selection of the ρ Oph Dark Cloud within the Upper Sco Region

In the main body of this paper, we select only stars within the actively star-forming region Ophiuchus. Here, we show the analysis results when selecting the entire Upper Sco region and discuss how the VSF of the stars depends on the stellar age.

The top panel of Figure A1 shows the spatial distribution of stars in Upper Sco and their age. Most of the stars outside Ophiuchus are older than 4 Myr, while stars inside Ophiuchus are much younger, primarily younger than 2 Myr. The SN explosion in this region traced by runaway stars took place ∼2 Myr ago (Neuhäuser et al. 2020), so the majority of stars in Upper Sco did not receive an energy injection from this event and thus would not exhibit characteristics of such energy injection in their VSF. Indeed, we see this trend in the lower panel of Figure A1. With the entire Upper Sco region taken into account, the VSF shows only a slight bump at ∼25 pc, while both the VSF of stars older than 5 Myr and of stars between 2 and 5 Myr do not show any discernible energy injection scale. Meanwhile, the VSF of stars younger than 2 Myr shows a clear energy injection scale at ∼25 pc, similar to what we obtain in the main paper. Because of these results, we elect to use only stars categorized as members of Ophiuchus, a majority of which are younger than 1 Myr, to analyze in the main body of this paper.

Appendix B: Effect of Galactic Rotation on the VSF of Hα

Figure B1 shows the VSF traced by Hα in four regions with and without Galactic rotation correction. In three out of the four regions, when a rotation profile is applied, the VSF does not significantly differ from the VSF without a rotation correction up to ∼20 pc. This behavior is expected because at small ℓ the rotation correction between two points is very similar. Thus, $\delta \vec{v}$ between two points essentially cancel out the rotation correction, and the shape of the VSF is preserved. For two points separated by a large distance, the rotation correction applied to them can differ more and hence results in a shift in the VSF, as Figure B1 shows. However, the overall shifts in the VSFs due to rotation correction are moderate and do not affect the conclusions of our work.

Appendix C: Effect of the Size and Location of the Hα Window

Another potential source of bias comes from the exact size and location of the window chosen to calculate the VSF of Hα within each region. In Section 3.2, we caution that the bumps seen in the VSF of Hα is a potential artifact because it is close to the size of the window. Such effect is noted as the apparent driving scale in simulations of multiphase turbulence and is equal to one-half of the box size (e.g., Mohapatra et al. 2022). To test whether this effect is also present in the VSF of Hα, we select a solid angle that encompasses all the stars in a region and compute the VSF within this window. We then compute the VSF for a window that is 4 deg² larger and repeat this a second time (see top panel in Figure C1).

**Figure C1.** Top panel: illustration of our Hα window tests. The red box is the smallest window within which all stars in Taurus belong. The two boxes in green (the window used in the main analysis) and light blue are each 4 deg² bigger than the former. Inside the red box, we also test the variation of the VSF within the region by dividing the window into four equal-sized boxes. Bottom-left panel: the VSF of Taurus traced by Hα, plotted across three different window sizes. Bottom-right panel: the Hα VSFs of Taurus in the four quadrants within the red window.
Download figure:
Standard image High-resolution image

In all four regions, the VSFs of Hα are independent of the window size up to ℓ ∼ L/2, where L is the window size. In Taurus and Perseus, the bump in the VSF shifts to larger ℓ when the window size increases. Thus, we conclude that in these regions, the VSF bump is likely caused by the box limit rather than a real physical driving scale. As an example, the bottom-left panel of Figure C1 shows the VSF of Hα in Taurus as a function of window size.

We also use Taurus as an example to illustrate the spatial variation of turbulence on small scales. We divide the smallest box (8 deg²) in the previous experiment into four equal-sized quadrants. The bottom-right panel of Figure C1 shows the Hα VSFs of these four quadrants in Taurus. Although the slopes of these four VSFs are roughly consistent with one another, their amplitudes vary by up to a factor of 3, indicating potentially large variations in turbulence between local patches of the ISM. The VSF of CO is in the best agreement with the bottom-left window (172° < l < 176° and −20° < b < −16°), where many CO cores were also identified (see Figure 2 of Qian et al. 2012 for detailed locations of CO cores in Taurus).

Turbulence in Milky Way Star-forming Regions Traced by Young Stars and Gas

Article metrics

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

Abstract

1. Introduction