Geometry, Kinematics, and Magnetization of Simulated Prestellar Cores

# Geometry, Kinematics, and Magnetization of Simulated Prestellar Cores

Che-Yu Chen11affiliation: Department of Astronomy, University of Virginia, Charlottesville, VA 22904, USA and Eve C. Ostriker22affiliation: Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA
###### Abstract

We utilize the more than 100 gravitationally-bound dense cores formed in our three-dimensional, turbulent MHD simulations reported in Chen & Ostriker (2015) to analyze structural, kinematic, and magnetic properties of prestellar cores. Our statistical results disagree with the classical theory of star formation, in which cores evolve to be oblate with magnetic field parallel to the minor axes. Instead, we find that cores are generally triaxial, although the core-scale magnetic field is still preferentially most parallel to the core’s minor axis and most perpendicular to the major axis. The internal and external magnetic field directions are correlated, but the direction of integrated core angular momentum is misaligned with the core’s magnetic field, consistent with recent observations. The ratio of rotational/total kinetic and rotational/gravitational energies are independent of core size and consistent in magnitude with observations. The specific angular momentum also follows the observed relationship L/M\propto R^{3/2}, indicating rotation is acquired from ambient turbulence. With typical E_{\mathrm{rot}}/E_{K}\sim 0.1, rotation is not the dominant motion when cores collapse.

magnetohydrodynamics (MHD) – stars: formation – turbulence

## 1. Introduction

In molecular clouds (hereafter MCs), multi-scale supersonic flows compress material and initiate creation of filamentary structures (André et al., 2014). Within filaments, some of the overdense regions will shrink under self-gravity to form prestellar cores and then collapse to create protostellar systems, which later become stars (Shu et al., 1987). Dense cores are therefore the immediate precursors of at least low-mass stars or close binary systems (McKee & Ostriker, 2007). Their properties provide the initial conditions of star formation, and determine the local environment of protostellar disks and outflows.

Cores are observed in dust continuum and molecular lines. Recent results from the Herschel Gould Belt Survey (André et al., 2010) suggest that dense cores are mostly associated with filaments (Könyves et al. 2010; or see review in André et al. 2014). This association is consistent with the theoretical expectation that thermally supercritical filaments (mass-per-unit-length M/L>2{c_{s}}^{2}/G; Ostriker 1964) would fragment longitudinally into cores (e.g. Inutsuka & Miyama, 1992, 1997). However, in a turbulent environment like a MC, dense filaments are not quiescent structures in which perturbations slowly grow. Rather, various simulations with turbulence have shown that filaments and cores develop simultaneously (e.g. H. Gong & Ostriker, 2011; Chen & Ostriker, 2014, 2015; Gómez & Vázquez-Semadeni, 2014; Van Loo et al., 2014; M. Gong & Ostriker, 2015), in contrast to the two-step scenario, because multi-scale growth is enabled by the non-linear perturbation generated by turbulence.

In combination with turbulence and gas gravity, magnetic effects are considered one of the key agents affecting the dynamics of star formation in MCs, at all physical scales and throughout different evolutionary stages (McKee & Ostriker, 2007). At earlier stages and on larger scales, the magnetic field can limit compression in turbulence-generated interstellar shocks that create dense clumps and filaments (Mestel & Spitzer, 1956). Meanwhile, the core-scale magnetic field is expected to be important in affecting the gas dynamics within cores, and is interconnected with the cloud-scale magnetic field. The magnetic field within collapsing cores provides the main channel for the gas to lose angular momentum via “magnetic braking” during the collapse of prestellar cores and the formation of protostellar disks (Mestel & Spitzer 1956; Gillis et al. 1974; Mouschovias 1976, 1991; see review in Li et al. 2014).

In strict ideal MHD, magnetic braking can be simply understood as the inner, faster-rotating material being slowed down by the outer, more slowly-rotating material because they are interconnected by magnetic field lines (Mouschovias & Paleologou, 1979, 1980). Numerical simulations also showed that the formation of rotationally supported disks is suppressed by catastrophic magnetic braking, unless the dense cores are weakly magnetized to an unrealistic level (Allen et al., 2003; Hennebelle & Fromang, 2008; Mellon & Li, 2008; Hennebelle et al., 2011). Many solutions have been proposed to solve this problem, including non-ideal MHD effects (Krasnopolsky et al., 2010; Li et al., 2011; Machida et al., 2011; Dapp et al., 2012; Tomida et al., 2013), turbulence-induced diffusion (Santos-Lima et al., 2012; Seifried et al., 2012, 2013; Joos et al., 2013), and the magnetic field-rotation misalignment (Hennebelle & Ciardi, 2009; Ciardi & Hennebelle, 2010; Joos et al., 2012; Krumholz et al., 2013). The initial magnetic field structure and strength within prestellar cores are therefore important for late evolution during core collapse, since they control the efficiency of this magnetic braking process.

Rotation is also important in the evolution leading to the creation of protostellar systems within dense cores. The angular momentum of star-forming cores is a critical parameter in protostellar evolution, but its origin is not well understood (see review in Li et al., 2014). It is known that some dense cores show a clear gradient in line-of-sight velocity, while others have a relatively random velocity field (e.g. Goodman et al., 1993; Caselli et al., 2002). The observed velocity gradient is commonly used, when present, to estimate a core’s angular momentum. Observational and theoretical understanding of core angular momentum is important as this property is essential to subsequent evolution: whether a single star or multiple system is formed (Tohline, 2002), and whether a large or small disk is produced (see review in Li et al., 2014).

There have been several observational projects aimed at resolving the velocity structure within dense cores. Linear fitting is generally applied to observed velocity gradients across cores, regardless of the complex nature of the velocity field. It is assumed that rigid body rotation applies and that the angular speed is roughly the gradient of line-of-sight velocity (Goodman et al., 1993). Previous observations (Goodman et al., 1993; Caselli et al., 2002; Pirogov et al., 2003; X. Chen et al., 2007; Tobin et al., 2011) have found a power-law relationship between the specific angular momentum, L/M, and radius, R, for dense cores/clumps with radii \sim 0.01-1 pc, L/M\propto R^{\alpha} with \alpha\approx 1.5.

The L/M-R correlation over a huge range of spatial scales suggests that gas motion in cores originates at scales much larger than the core size, or the observed rotation-like features may arise from sampling of turbulence at a range of scales (Burkert & Bodenheimer, 2000). Simulated cores from our previous work (Chen & Ostriker, 2014, 2015, hereafter CO14 and CO15) provide a suitable database to test whether the L/M-R correlation extends to even smaller scales (R \sim 0.008–0.05 pc), and whether the assumption of rigid body rotation is valid in simulated cores. These questions are one focus of the present study.

In this paper, we present our results on structural, kinematic, and magnetic properties of simulated prestellar cores formed in three-dimensional (3D) MHD simulations; the simulation ingredients include convergent flow, multi-scale turbulence, and gas self-gravity. This set of simulations was first introduced in CO15 to investigate the forming process of prestellar cores and how it correlates with the cloud environment. We have shown in CO15 that these cores have masses, sizes, and mass-to-magnetic flux ratios similar to the observed ones. Here, we extend our previous study to include analysis of the geometry and kinematic features of cores, as well as the relative direction of core-scale magnetic field.

The outline of this paper is as follows. We describe our simulation models, numerical methods, and analytic algorithms in Section 2. Our results on core geometry and other structural properties are presented in Section 3, while Section 4 focuses on the kinetic features of dense cores. Based on these results, we discuss the origin of core angular momentum in Section 5. Finally, Section 6 summarizes our conclusions.

## 2. Method

### 2.1. Numerical Simulations

The simulations we analyze here are described in CO14, CO15 and summarized here. These isothermal, self-gravitating ideal-MHD simulations consider the immediate surroundings of strongly-magnetized core-forming regions (1 pc{}^{3}) within MCs. For all simulations, the numerical resolution is 512^{3}, yielding cells of individual size \Delta x\approx 0.002 pc. The model establishes super-Alfvénic convergent flows (powered by the cloud-scale supersonic turbulence) compressing diffuse, turbulent gas to form denser, star-forming clumps in post-shock regions (see e.g. Figure 1). The idealized setup of a local converging turbulent flow gives us better control of the post-shock environment based on the simulation initial conditions as discussed in CO14, CO15, and Chen et al. (2017), which is crucial for analyzing the connection between the properties of prestellar cores and the environment they formed within.

The main model parameters (inflow Mach number {\cal M}_{0}\equiv v_{0,z}/c_{s} and background magnetic field strength B_{0}) are listed in Table 1. For all simulations, we initialize the background magnetic field in the x-z plane with an angle 20^{\circ} with respect to the converging flow along \pm\hat{z}. We consider two sets of parameters (M5, M10, M20 for varying {\cal M}_{0} and B5, B10, B20 for varying B_{0}, where M10 and B10 are the same model), which generates a range of core-forming environment and density structures in the post-shock region (Figure 1; also see Figure 3 of CO15). The post-shock conditions are determined by a combination of many mechanisms including the jump conditions of oblique MHD shocks (see CO14 and CO15) and the dynamics of secondary convergent flow (Chen et al., 2017).

Though we list the averaged post-shock plasma beta \beta_{\mathrm{ps}}\equiv 8\pi\rho_{\rm ps}{c_{s}}^{2}/{B_{\rm ps}}^{2} and Alfvén Mach number {\cal M}_{\mathrm{A,ps}}\equiv v_{\rm ps}/v_{\rm A,ps}=v_{\rm ps}\cdot\sqrt{4% \pi\rho_{\rm ps}}/B_{\rm ps} for each simulation model in Table 1, we note that one should not treat the post-shock environment using one single value of these parameters. This is especially true for those conditions that generate a dense sub-layer at the mid-plane of the post-shock region. Figure 1 shows the density (in \log cm{}^{-3}) and magnetic field directions (white streaklines) in slices cut through the center of each simulation box at y=0.5 pc. The large-scale converging flow is along the z axis. A dense, thin sub-layer can be clearly seen in models B20 and M5 but not in models B5 and M20, which generally agrees with the prediction in Chen et al. (2017). This cannot be derived from the averaged post-shock conditions (\beta_{\mathrm{ps}} and {\cal M}_{\mathrm{A,ps}}) listed in Table 1. As discussed in Chen et al. (2017), this sub-layer is created by supersonic secondary convergent flows and thus is relatively stagnant, which means dense clumps within it will undergo a more quiescent process to form prestellar cores. Roughly speaking, we expect the core-forming environment created in models B20 and M5 to be more dominated by the magnetic field, while gas turbulence plays a more critical role during core formation in models B5 and M20. We therefore rearrange the order of the models to group B20 and M5 as well as B5 and M20 in our following analysis and plots.

For each set of model parameters, we run 6 simulations with different realizations of the input turbulence. At the time the most evolved core collapses (when n_{\rm max}\geq 10^{7}~{}{\rm cm}^{-3}), we identify gravitationally bound cores as regions with E_{\rm grav}+E_{\rm thermal}+E_{\rm mag}<0 (see CO14, CO15 and below for details). This yields a total of 186 well-resolved self-gravitating cores within the post-shock layer.

### 2.2. Measuring Core Properties

For each simulation, we apply the GRID core-finding method to identify dense cores using the largest closed gravitational potential contours around single local potential minimums. The original GRID core-finding routine is developed and discussed in H. Gong & Ostriker (2011), while the MHD extension (to include measurement of cores’ mass-to-magnetic flux ratios) is implemented by and adopted in CO14 and CO15. Here, we further update the GRID core-finding method with angular momentum evaluation as well as derivation of rotational and turbulent energies. In this study, we also investigate the geometry of dense cores by applying the principal component analysis (PCA) to define the three axes of the core (a, b, and c from longest to shortest), and to calculate the aspect ratios (see e.g. H. Gong & Ostriker, 2011; M. Gong & Ostriker, 2015). Note that, similar to our previous studies (CO14, CO15), identified cores with less than 27 cells are not considered in the analysis. Figure 2 is a sample map showing the simulated gas structure (in column density integrated along z) and examples of identified cores from the GRID core-finding method.

### 2.3. Core Angular Momentum

For each core, it is straightforward to define and calculate the net angular momentum \mathbf{L} by integrating the relative angular momentum of each cell with respect to the center of mass over the whole volume:

 \mathbf{L}=\sum_{i}\rho_{i}\Delta V\cdot\left(\mathbf{r}_{i}-\mathbf{r}_{% \mathrm{CM}}\right)\times\mathbf{v}_{i} (1)

(\Delta V=\Delta x^{3} is the volume of a single simulation cell). This vector, together with the coordinate of the center of mass, determines the principal rotational axis for this core, \hat{\bf L}={\bf L}/L; here L=|{\bf L}| is the magnitude of the net angular momentum of the core. The total rotational inertia of the core around this rotational axis can be derived by firstly determining the projected radius for each cell:

 \mathbf{r}_{i,\perp}=\left(\mathbf{r}_{i}-\mathbf{r}_{\mathrm{CM}}\right)-% \left[\left(\mathbf{r}_{i}-\mathbf{r}_{\mathrm{CM}}\right)\cdot\hat{\mathbf{L}% }\right]\hat{\mathbf{L}} (2)

and then integrating over the whole volume:

 I\equiv\sum_{i}\rho_{i}\Delta V\cdot|\mathbf{r}_{i,\perp}|^{2} (3)

The mean angular velocity \Omega and the net rotational energy E_{\mathrm{rot}} of the core are

 \displaystyle\Omega \displaystyle\equiv L/I, (4) \displaystyle E_{\mathrm{rot}} \displaystyle\equiv\frac{1}{2}I\Omega^{2}. (5)

To compute the turbulent energy within cores under the assumption of rigid-body rotation, we must first calculate the velocity from rotation for each cell:

 \mathbf{v}_{\mathrm{rot},i}=r_{i,\perp}\Omega\left(\hat{\mathbf{L}}\times\hat{% \mathbf{r}}_{i,\perp}\right). (6)

The turbulent energy of the core under the assumption of rigid-body rotation is then

 E_{\mathrm{turb}}=\frac{1}{2}\sum_{i}\rho_{i}\Delta V\cdot\left(\mathbf{v}_{i}% -\mathbf{v}_{\mathrm{rot},i}\right)^{2}. (7)

The total kinetic energy is therefore defined as E_{\rm total}\equiv E_{\rm rot}+E_{\rm turb}.

### 2.4. The Ring-fit of Core Rotation

The simple assumption of rigid-body rotation adopted in Equations (6) and (7) is not always true. Here we consider a different approach of measuring the core’s rotational motion by binning the core into several “rings” and allowing individual rings to have different rotational speed.

The principal rotational axis given by \hat{\bf L} and \mathbf{r}_{\mathrm{CM}} defines a cylindrical coordinate system for the core: (\mathbf{r}_{i,\perp},h_{i}), where the “height” can be calculated as

 h_{i}=\left(\mathbf{r}_{i}-\mathbf{r}_{\mathrm{CM}}\right)\cdot\hat{\mathbf{L}}. (8)

We bin the cells by r_{i,\perp} and h_{i} to define local rings within the core; each ring is chosen to be 3-cell-wide in both radius and height. Similar to the cell-number requirement we applied to cores, if there are fewer than 27 cells assigned to a ring, that ring is not included in follow-up calculations. We also exclude cores with less than 4 rings (2 in both r and h directions) to improve the statistics.

For each ring (r,h), we calculate the net angular momentum (through its center of mass) \mathbf{L}_{\mathrm{ring}}(r,h), the rotational inertia I_{\mathrm{ring}}(r,h), the mean angular velocity \Omega_{\mathrm{ring}}(r,h), and the rotational energy E_{\mathrm{rot,ring}}(r,h) in the same way as Equations (1)-(5). The ring-fit mean angular velocity for a core is then defined as the inertia-weighted average among \Omega_{\mathrm{ring}}(r,h):

 \Omega_{\mathrm{ring}}=\frac{\sum\limits_{(r,h)}I_{\mathrm{ring}}(r,h)\cdot% \Omega_{\mathrm{ring}}(r,h)}{\sum\limits_{(r,h)}I_{\mathrm{ring}}(r,h)}. (9)

The total rotational energy derived from the ring fit is

 E_{\mathrm{rot,ring}}=\sum_{(r,h)}E_{\mathrm{rot,ring}}(r,h). (10)

where E_{\rm rot,ring}(r,h)=(1/2)I_{\rm ring}(r,h)\Omega_{\rm ring}(r,h)^{2}. The turbulent energy from the ring fit E_{\mathrm{turb,ring}} is obtained by following Equations (6) and (7), but with \Omega_{\mathrm{ring}}(\mathbf{r}_{i,\perp},h_{i}) instead of \Omega and \hat{\mathbf{L}}_{\mathrm{ring}}(\mathbf{r}_{i,\perp},h_{i}) instead of \hat{\mathbf{L}} in Equation (6). The total kinetic energy from the ring fit is E_{\rm total,ring}\equiv E_{\rm rot,ring}+E_{\rm turb,ring}. By comparing the ratios E_{\mathrm{turb}}/E_{\mathrm{total}} and E_{\mathrm{rot}}/E_{\mathrm{total}} to E_{\mathrm{turb,ring}}/E_{\mathrm{total,ring}} and E_{\mathrm{rot,ring}}/E_{\mathrm{total,ring}}, we can infer whether rigid-body rotation is a good approximation for our simulated cores (see Section 4.2).

Note that, because of the selection criteria, the number of cores for which we performed this ring-fit and energy analysis is less than the number of cores considered in the rigid-body rotation analysis discussed in Section 2.2 (see Tables 1 and 2). Also, the number of cells considered in each ring-fit might be less than the total number of cells within a core, which makes it inappropriate to directly compare the absolute values of E_{\mathrm{rot}} and E_{\mathrm{rot,ring}}, E_{\mathrm{turb}} and E_{\mathrm{turb,ring}}, or E_{\mathrm{total}} and E_{\mathrm{total,ring}}. We therefore use the relative values, E_{\mathrm{rot}}/E_{\mathrm{total}} and E_{\mathrm{rot,ring}}/E_{\mathrm{total,ring}} when cross-comparing these two fitting methods (see Table 2).

## 3. Core Geometry and Physical Properties

Table 1 summarizes the physical properties measured from simulated cores. Below we discuss the orientation (Section 3.1), geometry (Section 3.2) and the relative orientation of the magnetic field (Section 3.3) among these cores and between different models. Kinematic features related to the core’s angular momentum are discussed in Section 4. For reference, we note that (1) \hat{z} is the direction of converging flow and the “small” dimension for the post-shock layer within which; (2) \hat{x} is the primary direction of the magnetic field in the post-shock layer (see Figure 1); (3) the filaments within which cores form are in the \hat{x}-\hat{y} plane (see Figure 2 and CO15).

### 3.1. Core Orientation

For each core there are three PCA-defined axes/radii a, b, and c (from longest to shortest), and we can measure their directions by calculating the relative values of their x-, y-, and z-components, denoted as |a_{x}/a|, |a_{y}/a|, |a_{z}/a|, and so on. Figure 3 depicts the direction of the major and minor axes of all cores by showing the scatter plots of |a_{x}/a| vs. |a_{y}/a| (left) and |c_{x}/c| vs. |c_{z}/c| (middle); also shown is the scatter plot of |B_{x}/B_{\mathrm{tot}}| vs. |B_{z}/B_{\mathrm{tot}}|, the normalized x- and z-components of the average magnetic field within individual cores. We find that most cores have their major axes lying in the x-y plane of the simulation box ({|a_{x}/a|}^{2}+{|a_{y}/a|}^{2}\approx 1), with minor axes in the x-z plane ({|c_{x}/c|}^{2}+{|c_{z}/c|}^{2}\approx 1). Only cores formed in the more turbulent environment (models B5 and M20) can have relatively large values of |a_{z}/a| or |c_{y}/c| ({|a_{x}/a|}^{2}+{|a_{y}/a|}^{2}<1 or {|c_{x}/c|}^{2}+{|c_{z}/c|}^{2}<1). The result that the major axis lies in the x-y plane is not surprising as the converging flow along \hat{z} creates a dense post-shock layer, within which filaments and then cores form, in the x-y plane. Especially in the case of core formation in a stagnant sub-layer, it is reasonable to expect the major axis to be perpendicular to the inflow direction. Furthermore, since it is difficult to compress the magnetic field (which is roughly in the x-direction in the post-shock region) in the y-direction, minor axes are likely to lie in the x-z plane.

Similarly, the right panel of Figure 3 shows that core-scale magnetic field tends to lie in the x-z plane, with ({B_{x}}^{2}+{B_{z}}^{2})/{B_{\mathrm{tot}}}^{2}\approx 1. More importantly, there is a systematic variation of the magnetic field direction between different models, from preferably along z (|B_{z}/B_{\mathrm{tot}}|\approx 1) for models M5 and B20, to preferably along x (models B5 and M20). This result is consistent with the magnetic field direction in the post-shock region; for oblique MHD shocks propagating along z, the direction of post-shock magnetic field B^{\prime} has direction B^{\prime}_{x}/B^{\prime}_{z}=r_{B}B_{0,x}/B_{0,z} relative to the pre-shock magnetic field B_{0} (see Equation (7) in Chen et al., 2017), where r_{B} is the compression ratio of B_{x}. For strong shocks, r_{B}\sim v_{0}/v_{\mathrm{A}x0} (see derivations in CO14), which means B^{\prime}_{x}/B^{\prime}_{z}\propto v_{0}/(B_{0,x}B_{0,z})\propto v_{0}/{B_{0% }}^{2}. Models with weaker inflow Mach number (M5) or stronger pre-shock magnetic field (B20) therefore have post-shock magnetic field better aligned along z, while models with stronger inflow Mach number (M20) or weaker pre-shock magnetic field (B5) tend to have post-shock magnetic field almost parallel to x-direction. This suggests that, regardless the presence of the stagnant sub-layer, cores barely alter the magnetic field structure before they reach the collapse stage.

### 3.2. Aspect Ratio

The three PCA-defined axis lengths a, b, and c yield three distinct aspect ratios, which provide a measurement of how spherical a core is. The ratio c/b indicates whether a core is more prolate (c/b\sim 1), while the ratio b/a\sim 1 indicates an oblate core. Cores with c/b<1 and b/a<1 are triaxial. The median values of these ratios of each simulation model are listed in Table 1. Note that “cores” with b/a<0.2 (i.e. the major axis is at least 5 times longer than the other two axes) are considered elongated filaments and are not included in this study.

Figure 4 shows a scatter plot (c/a vs. b/a; top panel) and a histogram (c/b; bottom panel) of the ratios between the three axis lengths of simulated cores. Looking at the top panel of Figure 4, these cores have a wide range of axis ratios, and in general fall within the “triaxial” regime. Among cores from different models, models with stronger turbulence (M20) or weaker magnetization (B5) seem to have more prolate cores (closer to the diagonal b=c), while models with weaker turbulence (M5) or stronger magnetization (B20) tend to have low c/a (<0.5) and are more triaxial. This tendency can also be clearly observed in the bottom panel of Figure 4, where the median value of c/b shifts from \sim 0.5 for models M5 and B20, to \sim 0.6 for models M10B10 and M20, to \sim 0.7 for model B5.

We argue that this difference in aspect ratios of dense cores is caused by the existence of a stagnant sub-layer in some environments. As we have shown in Figure 3, cores formed within the stagnant sub-layers (models M5 and B20) tend to have their minor axes c along \hat{z}, because the stagnant sub-layers put limits on the core growth along the z direction by its thickness (see Figure 1). In this situation core formation is roughly two-dimensional in the x-y plane, within the sub-layer. On the other hand, in the situation without the presence of a stagnant sub-layer, local turbulence and the post-shock gas flow (which roughly follows the post-shock magnetic field direction on the x-z plane; see discussion in Chen et al. 2017) will lead to a more perturbed core-forming process even though the magnetic pressure is dominant in the post-shock region. We note that our simulated cores typically have b/a\sim 0.6 and c/a\sim 0.3 (see Table 1), inconsistent with the classical picture, in which cores are expected to be oblate (a\sim b\gg c) in strongly-magnetized regions (like the post-shock regions of our simulations).

Interestingly, the aspect ratio of dense cores does not seem to be strongly affected by the magnetization of the cores themselves. Figure 5 shows the scatter plots of the three aspect ratios b/a (top panels), c/a (middle panels), and c/b (bottom panels) as functions of core mass (left panels) and normalized mass-to-magnetic flux ratio \Gamma (right panels). The distributions of all three aspect ratios seem totally random with respect to \Gamma. This suggests that the magnetization level within cores is not a dominant factor controlling core geometry, at least at the evolutionary stage we measure. In contrast, there is a tendency for more massive cores to have smaller c/a and c/b, and for less massive cores to be more prolate (c\sim b). In detail, we find that while the overall core volume increases \propto M^{3/2} (due to the M\propto R^{2} correlation discussed in CO15), the range of a exceeds the range of c (and d\log a/d\log M>d\log c/d\log M at the high-mass end). Therefore, this may simply reflect the fact that the most massive cores have larger a.

### 3.3. Magnetic Field Orientation

It is interesting to investigate possible alignment between the magnetic field and core geometry. Figure 6 shows the histograms of the angle between the mean magnetic field within the core and its major (top panel) and minor (bottom panel) axes. These distributions show that cores have preferential alignments with respect to the local magnetic field, with minor axes preferentially more parallel to the field and major axes preferentially more perpendicular to the field. This is especially true in models with stronger pre-shock magnetic field (B20) or weaker pre-shock inflow Mach number (M5), which results in a less-perturbed post-shock environment, including the stagnant sub-layer, in which the cores form. However, we note that even though the longest axis is preferentially perpendicular to the magnetic field, cores are not strictly oblate (a\sim b), as would be expected for the very strongly magnetized case based on the classical picture of star formation. Also, there is no evidence that more oblate cores have short axes better aligned with the magnetic field, as would be expected in the classical picture. On the other hand, only in models without the stagnant sub-layers in the post-shock regions (highly turbulent (M20) or weakly magnetized (B5) models) do we find magnetic fields almost perpendicular to the core minor axes (large \measuredangle[B, c]), because the strong velocity turbulence can interfere the process of anisotropic condensation along magnetic field lines. However, we caution that cores tend to be more prolate (b\sim c) in these models (see Figure 4), and therefore the direction of c itself is less significant in this situation.

## 4. Kinematic Features

### 4.1. Integrated Angular Momentum and Rotational Energy

The median values of angular velocity \Omega and rotational energy E_{\mathrm{rot}} calculated from the integrated angular momentum L (see Section 2.2) are listed in Table 2. Also listed are the turbulent energy E_{\mathrm{turb}} derived using Equation 7, and the resulting rotational energy ratio E_{\mathrm{rot}}/E_{\mathrm{total}}. In all models, E_{\mathrm{rot}}/E_{\mathrm{turb}}\sim 0.1, and E_{\mathrm{rot}}/E_{\mathrm{total}}\sim 0.01-0.1. This suggests that rotation is not the dominant motion within prestellar cores. We discuss this further in Section 5.

#### 4.1.1 The Specific Angular Momentum and Core Geometry

As discussed in Section 2.2, the net angular momentum of a core, L, can be directly calculated by integrating through the core. When considering the magnitude of a dense core’s angular momentum, it is more common to adopt the specific angular momentum (the angular momentum per unit mass), L/M, than the net angular momentum itself (e.g. Goodman et al., 1993). Figure 7 shows the scatter plots of the aspect ratios b/a (top panel) and c/a (bottom panel) versus the specific angular momentum. In general, both aspect ratios decrease with increasing L/M; this can be interpreted as faster-rotating cores being more elongated. This is opposite to the naive expectation in which faster-rotating cores are more flattened/oblate (c/a\ll 1, b/a\sim 1), similar to the case for more rapidly rotating stars or planets. However, if we consider the definition of the net, integrated angular momentum (see Equation (1)), Figure 7 may simply reflect the fact that \Delta L_{i}\propto(r_{i}-r_{\mathrm{CM}}) (for given density and velocity). More specifically, for cores with similar mass, volume, and angular velocity, those with shapes more prolate will have larger angular momentum. This is one example of the potential risks of using the integrated angular momentum as the measurement of core rotation, because this value can be easily affected by the core geometry.

#### 4.1.2 Relative Orientation of the Angular Momentum

Similar to Figure 6, Figure 8 shows the histogram of the relative alignment of the integrated angular momentum \mathbf{L} with respect to the core’s major (top) and minor (bottom) axes. Though the core’s net angular momentum direction appears to be perpendicular to its major axis a, it seems to have no correlation with the minor axis c. This again is inconsistent with the naive expectation for a rotating oblate spheroid. We discuss this further in Section 5.

#### 4.1.3 Rotation-Magnetic Field Misalignment

In classical theory, the rotational axis is expected to be aligned with the magnetic axis of cores because magnetic braking is faster in perpendicular compared to parallel configurations (Mouschovias, 1979; Mouschovias & Paleologou, 1979). One of the most important breakthroughs in recent observations of magnetic field morphology within dense cores is that the magnetic field may not be aligned with the rotational axis as in classical theory (Hull et al., 2013, 2014). Our simulated cores provide the proper database for examining the rotation-magnetic field alignment in prestellar cores that are formed from a turbulent medium. Figure 9 illustrates the histogram (or the probability distribution function, PDF, top) and the cumulative distribution function (CDF; bottom) of the relative angles between \mathbf{B} and \mathbf{L} for all cores formed in our simulations. Though the PDF of each model show variations across [0^{\circ}, 90^{\circ}] and peak at different angles for different models, the CDF of all cores from various models combined appears to be a relatively straight line (random distribution) between [0^{\circ}, 90^{\circ}]. This agrees with the TADPOL result (Hull et al., 2013, 2014), which is shown in the bottom panel of Figure 9 as a step function, and suggests that there is no preferred orientation of core’s angular momentum with respect to its magnetic field. Though this picture differs from the classical model in which the rotational axis is aligned with the magnetic field, the rotation-magnetic field misalignment is critical in solving the magnetic braking catastrophe in protostellar disk formation (see review in Li et al., 2014).

The individual PDF of \measuredangle[B, L] for each model (Figure 9, top) shows that the relative angle between a core’s magnetic field and integrated angular momentum depends on the turbulence strength and magnetization level of the background cloud where the core forms. By looking at models with increasing pre-shock inflow Mach numbers (M5, M10, M20), we see that the median value of \measuredangle[B, L] shifts from small to large angles, which is also obvious for models with decreasing pre-shock magnetic field magnitude (B20, B10, B5). This is correlated with the turbulence level in the post-shock medium where the cores are formed, which is determined by the existence of the sub-layer; if the secondary convergent flow in the post-shock region is strong enough to create the sub-layer (as the cases of models B20 and M5), cores will form within this stagnant, dense sheet where magnetic field plays a more significant role than turbulence. This relatively quiescent core forming process is similar to the classic model, which results in nearly-aligned magnetic field and angular momentum within cores. On the other hand, in models B5 and M20 there is no sub-layer where the post-shock gas flow velocities (v^{\prime}\sim 0.1-0.2 km/s; see Equation (10) in Chen et al. 2017) collide and cancel before dense cores form, and therefore the post-shock gas momentum is applied directly to dense clumps as turbulence during the core forming process. Since the post-shock regions in these models are sub-Alfvénic, the gas flows near forming cores are likely along the magnetic field lines; if these anisotropic gas flows are the main sources of core’s angular momentum, it is not surprising that \measuredangle[B, L] is relatively large in these models.

Figure 10 compares the relative angles between L and a, c with the relative angles between B and L. Interestingly, we find that the distribution of \measuredangle[La] vs. \measuredangle[B, L] is concentrated in the upper-right half of the plot, meaning that the angles \measuredangle[La] and \measuredangle[B, L] cannot both be small (\lesssim 45^{\circ}). Correspondingly, the distribution of \measuredangle[Lc] vs. \measuredangle[B, L] appears to be a diagonal band from lower-left to upper-right of the plot, which means that these two angles might be positively correlated. These results agree with Figure 6, showing that the average magnetic fields within dense cores tend to align with the minor axes.

### 4.2. Non-rigid-body Rotation

Table 2 lists the average angular velocity \Omega_{\mathrm{ring}} and the relative rotational energy E_{\mathrm{rot,ring}}/E_{\mathrm{total,ring}} derived from our ring-fit method described in Section 2.4. Strikingly but not surprisingly, the rotational energy ratio from our ring-fit method is around a factor of 2 larger than that measured under the assumption of rigid-body rotation. Though this does not guarantee that the ring-fit method is more accurate for measuring core’s angular momentum (see the discussions in the Appendix), the fact that E_{\mathrm{rot}}/E_{\mathrm{total}} is only \sim 10~{}\% in either fitting method provides strong evidence that rotation is not the dominant motion within prestellar cores. We provide further discussion below in Section 5.

## 5. Discussion: The Origin of Core Angular Momentum

It has been known that the specific angular momentum of observed dense cores/clumps are correlated with their sizes, approximately following a power law, L/M\propto R^{\alpha}, with \alpha\approx 1.5. This was first reported in Goodman et al. (1993), and was later confirmed by many follow-up studies (e.g. Caselli et al., 2002; Pirogov et al., 2003; X. Chen et al., 2007). Also, these observations suggest that the ratio between rotational energy and gravitational energy, \beta_{E}\equiv(L^{2}/(2I))/(qGM^{2}/R) (where q=3/5 for a uniform density sphere; see definition in Goodman et al. (1993)), is relatively independent of core/clump size.

In Figure 11 we show the L/M vs. R relationship (bottom panel) from several observational studies (Goodman et al., 1993; Caselli et al., 2002; Pirogov et al., 2003; X. Chen et al., 2007; Tobin et al., 2011). We compare to the L/M vs. R relationship in our simulated cores considering R\equiv\sqrt[3]{a\cdot b\cdot c} (i.e. the geometric mean of the three axes), which follows the same correlation as the observations. The agreement of the simulations with observations confirms the positive correlation between the specific angular momentum and spatial scale of dense cores/clumps. Also plotted is the rotational-to-gravitational energy ratio \beta_{E} (middle panel) as a function of core radius R for both our simulations and observations. The independence of core rotational-to-gravitational energy ratio with core size is also confirmed. Quantitatively, the range of \beta_{E} for our simulated cores agrees with observations.

The fact that L/M\sim R\cdot v_{\mathrm{rot}}\propto R^{3/2} suggests that v_{\mathrm{rot}}\propto R^{1/2}. In combination with the well-known result that turbulent velocities increase roughly \propto R^{1/2} in supersonic turbulence (both in observations and simulations; see review in McKee & Ostriker 2007), this suggests that the rotational velocity in cores is inherited from the overall turbulent cascade. In addition, the top panel of Figure 11 shows the rotational-to-kinetic energy ratio, E_{\mathrm{rot}}/E_{K}, as a function of core radius R for both simulated and observed cores.111To estimate the kinetic energy of observed cores, we used the total observed linewidths \sigma_{v} within cores reported by the cited observation studies and calculated E_{K}\approx 3/2\cdot M_{\mathrm{core}}{\sigma_{v}}^{2}. For both simulated and observed cores, the range of E_{\mathrm{rot}}/E_{K} is similar and independent of R. This suggests that whatever size a core/clump is, it is sampling from the turbulence at that corresponding scale in setting its rotation. Though this rotational energy could be sub-dominant at core scales, it is essential to subsequent disk formation. Once a star-disk system forms in the interior (at much smaller scales than that studied here), the core’s velocity structure (and angular momentum) may be altered by outflows and/or other feedback mechanisms (see e.g. Offner & Chaban, 2017).

To investigate the accuracy of our estimated core angular momentum, we ran a set of tests to examine the analysis method, which is described in the Appendix. These tests show that the measured E_{\mathrm{rot}}/E_{K} ratio within the core could in principle reflect the relative significance of turbulence with respect to rigid-body rotation (see Figure 12, right panels). More importantly, we showed that for a pure-turbulent core (net angular velocity \Omega=0), the “projection” of turbulence within it will naturally lead to E_{\mathrm{rot}}/E_{K}\sim 0.1 (Figure 12, left panels), consistent with the values measured from our simulated cores. We therefore conclude that rotation is not the dominant motion within prestellar cores, and the ratio between the turbulence amplitude \sigma_{v} and maximum rotational speed (v_{\mathrm{rot,max}}\sim\Omega\cdot R_{\mathrm{core}}) must be \gtrsim 1 within prestellar cores.

## 6. Summary

In this paper, we investigated the >100 dense cores formed naturally in the CO15 MHD simulations to examine the structural, magnetic, and kinetic properties of prestellar cores with masses M_{\mathrm{core}}\sim 0.01-5~{}\mathrm{M_{\odot}} and sizes R_{\mathrm{core}}\sim 0.005-0.1 pc. We found that our simulated cores are generally triaxial, unlike the idealized oblate cores of classical theory. We showed that environmental effect plays an important role in shaping prestellar cores, especially by providing spatial constraints via ram or magnetic pressure. In addition, the formation of prestellar cores is strongly affected by gas turbulence, in the way that cores acquire rotational energy from local turbulence, which leads to the misalignment between magnetic fields and rotational axes within dense cores.

Our main conclusions are as follows:

1. When present, a stagnant sub-layer (see discussions in Chen et al., 2017) in the post-shock region (Figure 1) is critical in setting up the environment wherein prestellar cores form. Core formation within this sub-layer (models M5 and B20) is more quiescent and more similar to classical theory, while cores formed without this sub-layer (models M20 and B5) are more disturbed by local gas turbulence.

2. Cores preferentially have their major axes in the plane parallel to the shock front (x-y plane; see Figure 3, left), because the ram pressure of inflow limits core growth along z. This might help explain the mass-size relation of dense cores reported in both numerical (CO15) and observational (Kirk et al., 2013) studies, M\propto R^{k} with k\sim 2-2.5, because for cores formed within shock-compressed, locally-flat regions core growth is basically two-dimensional. On the other hand, we find that most of the cores have both their minor axes and mean magnetic fields lying in the x-z plane (Figure 3, middle and right panels) defined by the direction of inflow and the cloud-scale magnetic field. The minor axis is rarely along the \hat{y} direction because it is difficult for the gas to flow or compress the magnetic field in this direction.

3. Though cores are generally triaxial (Figure 4, top) rather than having the oblate shape often adopted in classical theory, the core-scale magnetic field is still generally aligned with core’s minor axis and perpendicular to the major axis (Figure 6). Only those cores formed under a more perturbed environment (without the stagnant sub-layer; models B5 and M20) can have magnetic field nearly perpendicular to their minor axes. However, these cores also tend to be more prolate (b\sim c; see Figure 4, bottom), which means the direction of c is less meaningful in these cases.

4. The integrated angular momentum vector within cores does not have a preferred orientation with respect to the core’s three axes (Figure 8), except being generally perpendicular to the major axis, which is a mathematical result from the definition of angular momentum (for cell i at a distance of r_{i} from the rotational axis, its angular momentum L_{i}\propto r_{i}). More importantly, there is no preferred alignment between the magnetic field and angular momentum within cores (Figure 9), as reported in the TADPOL observational survey (Hull et al., 2014) and a follow-up numerical study considering different viewing angles of two simulated protostellar envelopes (Lee et al., 2017). Since misalignment between a core’s rotational axis and magnetic field may be critical in reducing magnetic braking during core collapse, this may be important to understanding of protostellar disk formation.

5. Our analyses indicate that the commonly-adopted assumption of rigid-body rotation may underestimate the rotational motion in most dense cores. We presented a new method of calculating core angular momentum, the ring-fit (see Sections 2.4 and 4.2), which gives a factor of 2 higher measurement of rotational energy (see Table 2). Our results also suggest that the measured angular momentum within cores could simply be from the projection of ambient turbulence at core scale (see Figure 12) as previously suggested by Burkert & Bodenheimer (2000).

6. With our detailed analysis of core-scale kinematics, we have revisited the specific angular momentum-size correlation of dense cores/clumps reported in many observations (Figure 11, bottom). Our simulated cores fit with the observational results well, and extend to smaller spatial scales. The correlation of L/M with R over two orders of magnitude in spatial scale suggests that “rotation” within these cores/clumps shares the same origin with the velocity scaling consistent with larger-scale turbulence. We find that the rotational-to-gravitational energy ratio \beta_{E} (Figure 11, middle) and the relative rotational energy E_{\mathrm{rot}}/E_{K} ((Figure 11, top) have similar ranges and are independent of scale for both simulated and observed cores. Taken together, these results suggest that prestellar cores inherit their original angular momentum from cloud-scale turbulence, which may in part be driven by feedback (outflows, etc.) from other stars that formed earlier.

We note that our simulations and analyses only focus on cores in early (prestellar) evolutionary stages and therefore do not include effects from stellar feedback. There are also various idealizations in our simulations that could potentially affect our conclusions. The converging-flow setup intrinsically excludes scales \gtrsim L_{\mathrm{box}}, and therefore cannot capture effects of large-scale turbulence (including turbulence driven by feedback in other nearby stars) in development of core rotation. Also, we assumed uniform magnetic fields in the initial conditions, and this lack of magnetic field variation may affect the structure of local turbulence and hence rotation at core scales. Our results could also be biased by the limited parameter range (in terms of magnetic field strength, gas density, inflow velocity, etc.) investigated in this study. Nevertheless, the connection between core-scale angular momentum and the immediately-surrounding cloud-scale turbulence is clear in both our numerical results and previous observations.

We thank the referee for a very helpful report. C.-Y.C. is grateful for the support from Virginia Institute of Theoretical Astronomy (VITA) at the University of Virginia through the VITA Postdoctoral Prize Fellowship, and the support partly from NSF grant AST-1815784. The work of E.C.O. was supported by grant 510940 from the Simons Foundation.

## Appendix A The Significance of Integrated and Ring-fitted Angular Momenta

Here we describe the numerical tests that we conducted to examine our method of measuring angular momentum, for both the traditional rigid-body fit and our newly-developed ring fit. We constructed an ellipsoid with uniform density, \rho=1 in code units, and size a=20 and b=c=10 cells. This “core” is assigned a rigid-body angular speed \Omega=0, 1, or 10 in code units along a, plus a perturbed random velocity field with average amplitude \sigma_{v} ranging from 0.1 to 500 in code units; this leads to a range of rotational speed within the core (v_{\mathrm{rot,max}}=10 for \Omega=1 and 100 for \Omega=10). We then applied the same analysis for cores identified in our full simulations and measure the angular momentum, rotational energy, and turbulent energy of the core, and compared these values directly to imposed values.

### A.1. The Dependence of Measured Angular Momentum on Turbulent Level

The direct results are shown in the left panel of Figure 12. The \Omega=0 case gives the inherent numerical error in our measurement of E_{\mathrm{rot}}/E_{\mathrm{total}}, which is \sim 10~{}\% for both rigid-body and ring-fit methods in this case (we will discuss the dependence of this numerical error on cell resolution in Section A.2). This error value represents the minimum of E_{\mathrm{rot}}/E_{\mathrm{total}} ratio that can be correctly measured in this core; in both the \Omega=1 and \Omega=10 cases the measured E_{\mathrm{rot}}/E_{\mathrm{total}} values are truncated by the numerical error at large \sigma_{v}. Similarly, there is a numerical error for the E_{\mathrm{turb}}/E_{\mathrm{total}} ratio (\sim 10^{-2}) at small \sigma_{v}, but only for the ring-fit method (see discussion below).

The results from the two test cases \Omega=1 and \Omega=10 can be combined by plotting against the relative turbulent amplitude, which is defined as \sigma_{v}/v_{\mathrm{rot,max}} (right panel of Figure 12). We clearly see that both energy ratios change significantly when the turbulent level is around the same value as the maximum rotational speed (\sigma_{v}/v_{\mathrm{rot,max}}=1, indicated by the vertical lines in both plots). In addition, when the turbulent amplitude is about 5 times higher than the maximum rotational speed (\sigma_{v}/v_{\mathrm{rot,max}}\gtrsim 5), the measured E_{\mathrm{rot}}/E_{\mathrm{total}} ratios from both fitting methods can only serve as upper limits on the imposed values.

On the other hand, when the turbulent amplitude is only about 0.1 of the maximum rotational speed (\sigma_{v}/v_{\mathrm{rot,max}}\lesssim 0.1), the measured E_{\mathrm{turb}}/E_{\mathrm{total}} ratios from ring fit are truncated by the numerical error and again only serve as upper limits of the real values. In contrast, the rigid-body fit is not affected by the same truncating value, and has fairly accurate E_{\mathrm{turb}}/E_{\mathrm{total}} ratios even beyond \sigma_{v}/v_{\mathrm{rot,max}}\lesssim 0.01. This is likely because the velocity field within the core is constructed from rigid-body rotation, and therefore the ring-fit method will not have better performance than the rigid-body fit, especially when the turbulence amplitude is small (i.e. more comparable with pure rigid-body rotation).

With these results, we can in principle use the E_{\mathrm{rot}}/E_{\mathrm{total}} and E_{\mathrm{turb}}/E_{\mathrm{total}} values measured from our simulated prestellar cores to estimate \sigma_{v}/v_{\mathrm{rot,max}} in cores. Unfortunately, the range of E_{\mathrm{rot}}/E_{\mathrm{total}} measured in our simulations is mostly outside the zone where we can precisely estimate the turbulent amplitude (see the shaded regions in Figure 12). Nevertheless, we note that generally \sigma_{v}/v_{\mathrm{rot,max}}\gtrsim 1 inside prestellar cores (for both our simulated cores and for observed cores). This suggests that rotation is not the dominant motion within prestellar cores, and the measured angular momentum is indeed a projection of the turbulent velocity inside cores.

### A.2. The Dependence of Measured Angular Momentum on Cell Resolution

The numerical errors in E_{\mathrm{rot}}/E_{\mathrm{total}} and E_{\mathrm{turb}}/E_{\mathrm{total}} depend on the grid resolution, i.e. number of cells inside the core. The test case described in Section A.1 has \sim 8000 cells within the core, which is similar to some of the simulated cores considered in our main study. We repeated the same process on cores with more (\sim 10^{6}, a=50 and b=c=40 cells) and less (\sim 1000, a=10 and b=c=5 cells; see Figure 13) cells, and the truncating value of the E_{\mathrm{turb}}/E_{\mathrm{total}} ratio measured from the ring-fit method is \sim 10^{-3} and \sim 5\times 10^{-2}, respectively. Since cores formed in our simulations have numbers of cells ranging from \sim 30 to \sim 40000, the numerical errors for those cores will therefore be similar to that in the test cases. Even if the error in E_{\mathrm{turb}}/E_{\mathrm{total}} is slightly worse than \sim 10^{-2}, it will not affect the conclusions we draw from this test.

We also compared the ring-fit method with the rigid-body fit based on their performances under different resolution. From Figures 12 and 13, we found that the accuracy of ring-fit method highly depends on the number of cells within cores, even though the E_{\mathrm{rot}}/E_{\mathrm{total}} ratios measured from ring-fit method are always higher than that from rigid-body rotation. At higher resolution (e.g. Figure 12), ring-fitted rotational energy ratio E_{\mathrm{rot}}/E_{\mathrm{total}} follows the theoretical values better than the rigid-body fit; however, at lower resolution (Figure 13), the ring-fitted energy ratios deviate further away from the theoretical values (and the deviation happens at smaller turbulent amplitude) than the rigid-body fit. The numerical error in E_{\mathrm{rot}}/E_{\mathrm{total}} (the truncating value) from ring fit is also much larger than that in rigid-body fit for less-resolved core. Therefore, we conclude that the ring-fit is more applicable than the rigid-body fit only in larger (better resolved) cores.

## References

• Allen et al. (2003) Allen, A., Li, Z.-Y., & Shu, F. H. 2003, ApJ, 599, 363
• André et al. (2010) André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102
• André et al. (2014) André, P., Di Francesco, J., Ward-Thompson, D., et al. 2014, Protostars and Planets VI, 27
• Burkert & Bodenheimer (2000) Burkert, A., & Bodenheimer, P. 2000, ApJ, 543, 822
• Caselli et al. (2002) Caselli, P., Benson, P. J., Myers, P. C., & Tafalla, M. 2002, ApJ, 572, 238
• Chen & Ostriker (2012) Chen, C.-Y., & Ostriker, E. C. 2012, ApJ, 744, 124
• Chen & Ostriker (2014) Chen, C.-Y., & Ostriker, E. C. 2014, ApJ, 785, 69
• Chen & Ostriker (2015) Chen, C.-Y., & Ostriker, E. C. 2015, ApJ, 810, 126
• Chen et al. (2016) Chen, C.-Y., King, P. K., & Li, Z.-Y. 2016, ApJ, 829, 84
• Chen et al. (2017) Chen, C.-Y., Li, Z.-Y., King, P. K., & Fissel, L. M. 2017, ApJ, 847, 140
• X. Chen et al. (2007) Chen, X., Launhardt, R., & Henning, T. 2007, ApJ, 669, 1058
• Ciardi & Hennebelle (2010) Ciardi, A., & Hennebelle, P. 2010, MNRAS, 409, L39
• Dapp et al. (2012) Dapp, W. B., et al. 2012, A&A, 541, A35
• Fissel et al. (2016) Fissel, L., Ade, P., Angile, F., et al. 2016, ApJ, 824, 134
• Gillis et al. (1974) Gillis, J., Mestel, L., & Paris, R. B. 1974, Ap&SS, 27, 167
• Gómez & Vázquez-Semadeni (2014) Gómez, G. C., & Vázquez-Semadeni, E. 2014, ApJ, 791, 124
• H. Gong & Ostriker (2011) Gong, H., & Ostriker, E. C. 2011, ApJ, 729, 120
• M. Gong & Ostriker (2015) Gong, M., & Ostriker, E. C. 2015, ApJ, 806, 31
• Goodman et al. (1993) Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528
• Hennebelle & Fromang (2008) Hennebelle, P., & Fromang, S. 2008, A&A, 477, 9
• Hennebelle et al. (2011) Hennebelle, P., et al. 2011, A&A, 528, A72
• Hennebelle & Ciardi (2009) Hennebelle, P., & Ciardi, A. 2009, A&A, 506, L29
• Hull et al. (2017) Hull, C. L. H., Girart, J. M., Tychoniec, Ł., et al. 2017, ApJ, 847, 92
• Hull et al. (2013) Hull, C. L. H., Plambeck, R. L., Bolatto, A. D., et al. 2013, ApJ, 768, 159
• Hull et al. (2014) Hull, C. L. H., Plambeck, R. L., Kwon, W., et al. 2014, ApJS, 213, 13
• Inutsuka & Miyama (1992) Inutsuka, S.-I., & Miyama, S. M. 1992, ApJ, 388, 392
• Inutsuka & Miyama (1997) Inutsuka, S.-i., & Miyama, S. M. 1997, ApJ, 480, 681
• Joos et al. (2012) Joos, M., Hennebelle, P., & Ciardi, A. 2012, A&A, 543, A128
• Joos et al. (2013) Joos, M., et al. 2013, A&A, 554, A17
• Kirk et al. (2013) Kirk, J. M., Ward-Thompson, D., Palmeirim, P., et al. 2013, MNRAS, 432, 1424
• Könyves et al. (2010) Könyves, V., André, P., Men’shchikov, A., et al. 2010, A&A, 518, L106
• Krasnopolsky et al. (2010) Krasnopolsky, R., et al. 2010, ApJ, 716, 1541
• Krumholz et al. (2013) Krumholz, M. R., et al. 2013, ApJ, 767, L11
• Lee et al. (2017) Lee, J. W. Y., Hull, C. L. H., & Offner, S. S. R. 2017, ApJ, 834, 201
• Li et al. (2011) Li, Z.-Y., et al. 2011, ApJ, 738, 180
• Li et al. (2014) Li, Z.-Y., Banerjee, R., Pudritz, R. E., et al. 2014, Protostars and Planets VI, 173
• Machida et al. (2011) Machida, M. N., et al. 2011, PASJ, 63, 555
• McKee & Ostriker (2007) McKee, C. F., & Ostriker, E. C. 2007, ARA&A, 45, 565
• Mellon & Li (2008) Mellon, R. R., & Li, Z.-Y. 2008, ApJ, 681, 1356
• Mestel & Spitzer (1956) Mestel, L., & Spitzer, L., Jr. 1956, MNRAS, 116, 503
• Mouschovias (1976) Mouschovias, T. C. 1976, ApJ, 207, 141
• Mouschovias (1979) Mouschovias, T. C. 1979, ApJ, 228, 159
• Mouschovias (1991) Mouschovias, T. C. 1991, ApJ, 373, 169
• Mouschovias & Paleologou (1979) Mouschovias, T. C., & Paleologou, E. V. 1979, ApJ, 230, 204
• Mouschovias & Paleologou (1980) Mouschovias, T. C., & Paleologou, E. V. 1980, ApJ, 237, 877
• Offner & Chaban (2017) Offner, S. S. R., & Chaban, J. 2017, ApJ, 847, 104
• Ostriker (1964) Ostriker, J. 1964, ApJ, 140, 1056
• Pirogov et al. (2003) Pirogov, L., Zinchenko, I., Caselli, P., Johansson, L. E. B., & Myers, P. C. 2003, A&A, 405, 639
• Planck Collaboration Int. XIX (2015) Planck Collaboration Int. XIX. 2015, A&A, 576, A104
• Santos-Lima et al. (2012) Santos-Lima, R., et al. 2012, ApJ, 747, 21
• Seifried et al. (2012) Seifried, D., et al. 2012, MNRAS, 423, L40
• Seifried et al. (2013) Seifried, D., et al. 2013, MNRAS, 432, 3320
• Shu et al. (1987) Shu, F. H., Adams, F. C., & Lizano, S. 1987, ARA&A, 25, 23
• Stone et al. (2008) Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
• Tobin et al. (2011) Tobin, J. J., Hartmann, L., Chiang, H.-F., et al. 2011, ApJ, 740, 45
• Tohline (2002) Tohline, J. E. 2002, ARA&A, 40, 349
• Tomida et al. (2013) Tomida, K., et al. 2013, ApJ, 763, 6
• Van Loo et al. (2014) Van Loo, S., Keto, E., & Zhang, Q. 2014, ApJ, 789, 37
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters