CLUSTERED STAR FORMATION IN MAGNETIC CLOUDS: PROPERTIES OF DENSE CORES FORMED IN OUTFLOW-DRIVEN TURBULENCE

Following Nakamura & Li (2008), we identify dense cores from the density data cube using a variant of the CLUMPFIND algorithm of Williams et al. (1994). We adopt a threshold density ρ_{th, min} = 4ρ₀ (∼10⁵ cm⁻³) and a maximum density ρ_{th, max} = 400ρ₀ (∼10⁷ cm⁻³, above which a Lagrangian particle is created). Thus, the cores identified from our simulation data correspond roughly to the dense cores observed in the dust continuum emission and high density molecular tracers such as N₂H⁺ (J = 1–0) and H¹³CO⁺ (J = 1–0). The density distribution between ρ_{th, min} and ρ_{th, max} is divided into 10 bins equally spaced logarithmically. We include only those spatially resolved cores containing more than 50 cells in the analysis to ensure that their properties are determined with reasonable accuracy. We also choose the cores that do not include particles. In this sense, our identified cores correspond to starless cores. The minimum core mass and radius identified through this procedure are M_{c, min} ≡ ρ_{th, min} × 50ΔxΔyΔz ∼ 0.05 M_☉ and ∼0.01 pc, respectively. We note that our core identification procedure is not exactly the same as those used in observational studies, which are typically based on position–position–velocity data cubes of molecular line emission or column density maps from dust continuum emission.

Before discussing the physical properties of the identified cores, we show how the magnetic field affects the global density distribution. Figure 1 compares the column density distributions along the y-axis for the three models with different initial magnetic field strengths, at a stage when the star formation efficiency (SFE) has reached 16%. Since the star formation is retarded by the magnetic field, the evolution time tends to be longer at a given SFE for a stronger magnetic field. For comparison, the positions of the formed stars and the identified dense cores projected on the x–z plane are indicated by the dots and triangles, respectively, in Figure 1. The initial magnetic field direction is parallel to the x-axis (the abscissa).

Zoom In Zoom Out Reset image size

Figure 1 indicates that the global density distribution depends on the initial magnetic field strength. For the non-magnetic model (model N1), the global density distribution shows a large-scale filamentary structure that contains many small fragments (Figures 1(a) and (d)). The formation of the large-scale filament is probably because the turbulence energy is initially highest on the largest scale and the large-scale compression may be further amplified by self-gravity to create the filamentary structure.

For the weakly magnetized case (model W1), the global density distribution does not clearly show a large-scale filamentary structure by the stage shown in Figures 1(b) and (e). The dense part appears to contain many smaller filaments or elongated fragments that are distributed almost independently of the initial magnetic field direction. This is different from the non-magnetic case. The difference may be because the local turbulent motions amplify random magnetic field components (see Section 4.2), which slow down the formation of the large-scale filamentary structure.

In contrast, in the presence of a moderately strong magnetic field (model S1), the cloud material condenses preferentially along the magnetic field lines into a large-scale filamentary structure that is nearly perpendicular to the initial magnetic field direction. Several less-dense filaments are also associated with the large-scale filament (hereafter the main filament). The less-dense filaments appear to be more or less along the initial magnetic field direction and to merge toward the dense parts of the main filament. The dense cores are distributed primarily along the main filament. The column density distribution is reminiscent of the hub-filament structures of star-forming regions pointed out by Myers (2009). The observed hub-filament structures may indicate the presence of a moderately strong, spatially ordered magnetic field, instead of the weak magnetic field tangled by supersonic turbulence. We will revisit this issue in Section 4.3, where we compare dust polarization maps derived from the simulation data for models W1 and S1.

As shown in Li & Nakamura (2006) and Nakamura & Li (2007), our model clump has reached a quasi-equilibrium state within one global free-fall time. After that, the statistical properties of the cores remain largely unaltered. Therefore, to compare the statistical properties of the identified cores, we choose the same stage as that of Figure 1, i.e., the stage when 16% of the total mass has been converted into stars. In this subsection, we compare the core properties for the three representative models, models N1, W1, and S1. The total numbers of identified cores are 160, 105, and 105 cores for models N1, W1, and S1, respectively.

3.3.1. Radius, Mass, and Velocity Dispersion

We now focus on the physical properties of the cores identified in the three representative models. To facilitate comparison with observations, we present the core properties in dimensional units. Figures 2(a)–(i) show the histograms of the core radius, mass, and velocity dispersion of the nonthermal component. The minimum, maximum, mean, and median values are summarized in Table 2. The core radius, R_c, is defined as the radius of a sphere having the same volume as the core. The distributions of core radius are broadly similar in all cases. They range from ∼0.01 pc to ∼0.05 pc, all peaking around 0.02 pc. The maximum and mean core radii are somewhat larger in model S1 (the moderately strong field case) than in the other two. The core mass distribution shows a similar trend. In all models, the core mass ranges from ∼0.07 M_☉ to ∼4 M_☉. The mean value is somewhat larger for model S1 than for the other two. The mean core mass and radius are larger in model S1 probably because the stronger magnetic field can support more mass against gravitational collapse.

Zoom In Zoom Out Reset image size

Table 2. Summary of the Physical Properties of the Identified Cores

Property	Minimum	Maximum	Meana	Median
Model N1
Rc (pc)	0.013	0.041	0.020 ± 0.005	0.020
Mc (M☉)	0.075	3.51	0.67 ± 0.67	0.40
ΔVNT (km s−1)	0.13	4.53	1.06 ± 0.77	0.86
Model W1
Rc (pc)	0.013	0.034	0.020 ± 0.006	0.019
Mc (M☉)	0.076	3.90	0.62 ± 0.77	0.31
ΔVNT (km s−1)	0.13	14.4	0.99 ± 1.66	0.47
Model S1
Rc (pc)	0.013	0.048	0.024 ± 0.008	0.023
Mc (M☉)	0.073	4.39	0.86 ± 1.03	0.47
ΔVNT (km s−1)	0.086	2.64	0.46 ± 0.52	0.31

Notes. The gas temperature is assumed to be T = 20 K, corresponding to the sound speed of 0.266 km s⁻¹. ^aWith standard deviation.

Download table as:  ASCIITypeset image

The magnetic effect is even more prominent in the distribution of core velocity dispersion. For the non-magnetic model (model N1), the velocity dispersion ranges from ∼0.1 km s⁻¹ to 4 km s⁻¹, peaking around 1 km s⁻¹ (corresponding to a Mach number of ∼4). Most of the cores have supersonic internal motions in this case, which appears inconsistent with the predominantly subsonic or at most transonic internal motions inferred observationally (e.g., André et al. 2007; Walsh et al. 2007; Maruta et al. 2010; see also Section 4.1). Even a relatively weak magnetic field appears to affect the internal motions inside the cores significantly. For model W1, the median velocity dispersion is estimated at ∼0.5 km s⁻¹, corresponding to a Mach number of about 2, although the mean is almost the same as that of model N1. Most of the cores in this case still have supersonic velocity dispersions.

In contrast, for model S1, the majority of the cores have a velocity dispersion close to, or less than, ∼0.3 km s⁻¹, which is comparable to the sound speed (c_s = 0.266 km s⁻¹). The velocity dispersion ranges from ∼0.05 km s⁻¹ to 2 km s⁻¹, peaking around the median of 0.3 km s⁻¹. In other words, the strong magnetic field tends to reduce the core internal motions significantly. The reduction is probably because the magnetic field cushions the converging flows in the supersonic turbulence that are largely responsible for the core formation. The magnetic effect is also evident in the virial state of the cores, as we show below.

3.3.2. Velocity-dispersion–Radius Relation

The internal motions of the observed dense cores are often measured using the velocity dispersions derived from molecular line data. Here, we compare the velocity dispersions for the three representative models with different field strengths to study the effects of the magnetic field.

Figure 3 plots the nonthermal components of three-dimensional velocity dispersions of the identified cores against the core radii for (a) model N1, (b) model W1, and (c) model S1. The best-fit power laws are given by (ΔV_NT/km s⁻¹) = (2.2 ± 1.9)(R_c/pc)^{0.19 ± 0.22} (${\cal R} = 0.023$), (ΔV_NT/km s⁻¹) = (3.8 ± 9.0)(R_c/pc)^{0.35 ± 0.61} (${\cal R} = 0.070$), and (ΔV_NT/km s⁻¹) = (1.1 ± 0.48)(R_c/pc)^{−0.34 ± 0.12} (${\cal R} = 0.29$) for models N1, W1, and S1, respectively, where ${\cal R}$ denotes the correlation coefficient. It shows that there is no clear correlation between the velocity dispersion and core radius. Our cores do not follow the power-law linewidth–radius relation derived by Larson (1981). The same conclusion holds for dense cores formed in even more strongly magnetized, magnetically subcritical clouds (see Figure 18 of Nakamura & Li 2008). The reason is probably that the dense cores created out of turbulent molecular gas are generally not in virial equilibrium (the equilibrium is implied by Larson's laws) and that their formation is strongly influenced by external turbulent flows (see below).

Zoom In Zoom Out Reset image size

Quantitatively, the velocity dispersions of the cores depend on the initial magnetic field strength, as shown already in Figures 2(g)–2(i). In the absence of magnetic field (model N1), about 40% of the cores have velocity dispersions larger than 1 km s⁻¹; the majority of the cores have supersonic internal motions. As mentioned above, even a relatively weak magnetic field can influence the internal motions significantly. In model W1, although many cores still have large velocity dispersions, the median value is less than that in model N1 by a factor of ∼2. In the stronger field case (model S1), the median and mean of the velocity dispersions are estimated at 0.31 and 0.46 km s⁻¹. About 40% of the cores have velocity dispersions smaller than the sound speed of 0.266 km s⁻¹; the majority of the cores have subsonic or at most transonic internal motions, which is more consistent with observations. We conclude that the magnetic field plays a significant role in reducing the internal motions of the cores.

3.3.3. Virial Parameter

The virial parameter, the ratio of the virial mass to core mass, is often used as a measure of the gravitational boundedness of dense molecular cloud cores, particularly in observational studies. In Figure 4, we plot the virial parameters against the core mass for the three models. Here, we calculate the virial parameter as

$$\begin{equation} \alpha _{\rm vir}= \frac{5\Delta V_{\rm 1D}^2 R_{\rm c}}{GM_{\rm c}} \, \end{equation} \tag{ 1 }$$

where ΔV_1D is the one-dimensional FWHM velocity width including thermal contribution and the virial mass of each core is obtained by assuming that it is a uniform sphere. The virial mass depends on the density distribution. For a centrally condensed sphere with ρ∝r⁻², the virial parameter is reduced by a factor of 5/3. Note that the effects of the nonspherical mass distribution appear to be small as long as the aspect ratios of the cores are not far from unity (Bertoldi & McKee 1992).

Zoom In Zoom Out Reset image size

Figure 4 indicates that the relationship between the virial parameter and core mass depends on the initial magnetic field strength. In the absence of a magnetic field (model N1), the virial parameters are typically very large and have a large scatter. The scatter in the virial-parameter–mass plot comes mainly from the large scatter in the velocity-dispersion–radius relation (see Figure 3). The virial parameters often exceed 10–10² even for relatively massive cores, indicating that, for most of the cores, the gravitational energy is much smaller than the internal kinetic energy. The virial parameters tend to increase with decreasing core mass, and the lower bound of the virial parameters follows roughly a power-law relation of α_vir∝M^−2/3_c.

For the model with a weak magnetic field (model W1), the virial-parameter–mass relation is qualitatively similar to that of the non-magnetized model (model N1). Although the fraction of cores with large virial parameters is somewhat small for the weak magnetic field case, the scatter in the virial-parameter–mass relation remains large. About 30% of the cores still have virial parameters larger than 10.

On the other hand, for the model with a moderately strong magnetic field (model S1), the identified cores tend to follow a single power-law relation of α_vir∝M^−2/3_c, with a rather small scatter. Furthermore, most of the cores (about 9% of the cores) have virial parameters smaller than 10. This trend is consistent with the initially magnetically subcritical case (see Figure 17 of Nakamura & Li 2008). It is also in good agreement with the scaling found by Bertoldi & McKee (1992) for nearby GMCs. According to Bertoldi & McKee (1992), an object confined by the ambient pressure has a virial ratio that is proportional to (M_c/M_J)^−2/3, where M_J is the Jeans mass defined by Equation (2.13) of Bertoldi & McKee (1992). The fact that our cores follow the same relation implies that the ambient pressure plays an important role in confining our cores as well. In the next subsection, we perform a detailed virial analysis of the cores to quantitatively clarify the role of the surface pressure term in core dynamics.

3.3.4. Virial Analysis

The virial theorem is useful for analyzing the core dynamical properties. Here, we perform a detailed virial analysis using the virial equation in Eulerian coordinates (e.g., McKee & Zweibel 1992; Tilley & Pudritz 2007; Dib et al. 2007):

$$\begin{equation} S=\frac{1}{2}\ddot{I} + \frac{1}{2}\frac{d}{dt}\int _S \rho r^2 \mbox{\boldmath $v$} \cdot \mbox{\boldmath $dS$} = U+K+M+W+S+F \, \end{equation} \tag{ 2 }$$

where I = ∫_Vρr²dV, U = 3∫_VPdV, K = ∫_Vρv²dV, $W=-\int _V \rho \mbox{\boldmath $r$} \cdot \nabla \Psi dV$, M = (1/8π)∫_VB²dV, $S= - \int _S [P\mbox{\boldmath $r$}+\mbox{\boldmath $r$}\cdot (\rho \mbox{\boldmath $r$}\mbox{\boldmath $r$})] \cdot \mbox{\boldmath $dS$}$, and $F=\int _S \mbox{\boldmath $r$}\cdot \mbox{\boldmath $T_M$} \cdot \mbox{\boldmath $dS$}$. The quantities I, U, K, W, S, M, F, and $\mbox{\boldmath $T_M$}$ denote, respectively, the moment of inertia, internal thermal energy, internal kinetic energy, gravitational energy, the sum of the thermal surface pressure and dynamical surface pressure, internal magnetic energy, the magnetic surface pressure, and the Maxwell stress-energy tensor. The internal thermal, kinetic, and magnetic energies are always positive. Other terms can be either positive or negative. In particular, the gravitational term can be positive in crowded environments where the background gravitational field dominates the core dynamics (see e.g., Ballesteros-Paredes et al. 2009). In fact, as shown below, several cores have positive gravitational terms in our simulations, although no cores have positive gravitational terms for the initially magnetically subcritical case where the turbulent motions are not that strong because of the effects of magnetic cushion (see Figure 16 of Nakamura & Li 2008). The second term on the left-hand side of Equation (2) denotes the time derivative of the flux of moment of inertia through the core boundary. As is the standard practice (Tilley & Pudritz 2007; Dib et al. 2007; Nakamura & Li 2008), we ignore the left-hand side of Equation (2) in our discussion and consider a core to be in virial equilibrium if the sum U + K + W + S + M + F = 0.

The equilibrium line is shown in Figure 5, where the sum of the surface terms (S+F) is plotted against the gravitational term (W), both normalized to the sum of the internal terms (U + K + M). For the surface term, the external kinetic term is generally much larger than the magnetic term (F).

Zoom In Zoom Out Reset image size

For all three models, the majority of the cores lie below the virial equilibrium line, and thus they are expected to be unstable to contraction. In addition, the surface term appears more important than the gravitational term for most of the cores. In the absence of the magnetic field, about 10%–20% of the cores have the normalized surface term of |S + F|/|U + K + M| ≳ 5, whereas no cores have the normalized surface term of |S + F|/|U + K + M| ≳ 5 in the presence of moderately strong magnetic field. The magnetic field tends to reduce the contribution of the surface terms to the core dynamical state, although the surface term is still dominant even in the presence of strong magnetic field.

The core virial state is in contrast to the initially magnetically subcritical case. For the initially magnetically subcritical case, for most of the cores, the surface term (S+F) and gravitational term (W) tend to be smaller than the sum of the internal terms, U + K + M (Nakamura & Li 2008), and thus most of the cores are distributed around the virial equilibrium line. In addition, more massive cores tend to be distributed below the virial equilibrium line, and the importance of self-gravity tends to depend on the core mass. The gravitational term tends to be more important for more massive cores for the initially magnetically subcritical case, implying that the gravitational collapse induced by ambipolar diffusion regulates the core evolution. In contrast, for the magnetically supercritical cases studied in this paper, most of the cores are distributed farther from the virial equilibrium line. The surface term (S+F) tends to be larger than the sum of the internal terms, U + K + M for most of the cores. The core virial state shows no clear dependence on the mass. The surface term, which is typically dominated by the ambient supersonic turbulence, is still important even for massive cores.

The above result is different from that of Tilley & Pudritz (2007). They found, in their ideal MHD simulations, a larger contribution from the surface term to the virial equation for a stronger magnetic field, which is opposed to the trend we found. This difference likely comes from the fact that they only followed the very early phase of cloud evolution, when the initial supersonic turbulence has decayed significantly but star formation has yet to set in. Most of our cores are formed at much later times, after the decayed turbulence has been replenished by outflow feedback. Given this difference, it is perhaps not too surprising that the surface terms are much larger in our models than in theirs. This is particularly true for the weaker magnetic field case, because of a more rapid and violent star formation.

In our simulations, most of the identified cores have densities smaller than the critical density beyond which the Jeans condition is violated. For example, the Jeans condition is satisfied for 89% (model N1), 85% (model W1), and 95% (model S1) of the cores presented in this subsection. For the stronger initial magnetic field, the mean density of the cores tend to be smaller due to the magnetic support, and therefore the Jeans condition is satisfied for a larger number of cores. We also performed a run with 512³ for model S1, with the threshold density fixed, so that the Jeans condition is always satisfied in the entire computation box. To compare directly with the results presented in this subsection, we first resized the data obtained from the 512³ run to the 256³ grid data, and then applied the clumpfind method to the regridded data. We confirmed that the statistical properties of the identified cores are essentially the same as those presented in this subsection. For example, for all the cores identified from the 512³ data, the nonthermal three-dimensional velocity dispersions of the identified cores stay below 1 km s⁻¹, and its mean is estimated to be 0.43 km s⁻¹, comparable to that of the 256³ data, 0.46 km s⁻¹.

After exploring the effects of the magnetic field on the core properties, we now turn to the outflow strength. In our model, the strength of the outflow feedback is specified by the dimensionless parameter, f. This parameter is highly uncertain because it is generally difficult to determine the physical properties of molecular outflows accurately from observational data. In this paper, we adopt f = 0.4 as a fiducial value, following Matzner & McKee (2000). Recently, Nakamura et al. (2011b) discovered many molecular outflow lobes toward an extremely young cluster-forming clump, Serpens South, on the basis of CO (J = 3–2) observations, and derived the physical properties of the identified molecular outflows. They roughly estimated f ∼ 0.1–0.3 for these outflows, assuming that the outflow gas is optically thin. If a substantial amount of the outflow gas is optically thick, the outflow parameter f would be larger. In the following, we compare the core properties of three models with different outflow strengths: f = 0, 0.1, and 0.4. The fiducial value for plasma β of 0.2 is adopted for these models.

In Figures 6(a)–(c), we compare the velocity-dispersion– radius relations for the moderately strong magnetic models with no outflow feedback (f = 0), weak outflow (f = 0.1), and strong outflow (f = 0.4). In the model with no outflow feedback (model S3), the number of identified cores is smaller than in the other two models, and the mean core mass tends to be slightly larger. This difference may be due to the fact that the ambient turbulent motions are weaker without any outflow feedback and thus the turbulent fragmentation does not form smaller-scale structures as efficiently. Figure 6 indicates that the velocity-dispersion–radius relation does not depend strongly on the strength of the outflow feedback. For all three models, there is no clear relationship between the velocity dispersion and the core radius. Thus, it appears difficult to discriminate the effects of protostellar outflows based on the velocity-dispersion–radius relation alone.

Zoom In Zoom Out Reset image size

The surface term-gravitational term plots in Figures 7(a)–(c) indicate that the dynamical states of the cores depend somewhat on the outflow strength. In the absence of the outflow feedback, more cores are located near the dashed line (where the two terms are equal). The gravitational term appears more important for this model than for the other two, although the surface term is clearly still important. This is due, at least in part, to the fact that the global gravitational infall can generate fast motions even in the absence of protostellar outflow-driven turbulence.

Zoom In Zoom Out Reset image size

It is difficult, however, to gauge the effects of the outflow feedback from the virial parameters of the cores, which are often used in observational studies of the molecular cloud cores in star-forming molecular clouds. The difficulty is illustrated in Figures 8(a)–(c), where we show the virial-parameter–mass relations for the three models with different outflow strengths. In all cases, the distributions of the virial parameters are similar, following a power law of α_vir∝M^−2/3_c with rather small dispersions.

Zoom In Zoom Out Reset image size

In this subsection, we investigate how well the outflow axes are aligned with the global or initial magnetic field direction in parsec-scale cluster-forming clumps. In our simulations, the outflow direction is set to the local magnetic field direction at the position of each formed star. This should be a reasonable approximation since the outflows are most likely driven along the local magnetic field by field line twisting due to disk rotation (e.g., Matsumoto et al. 2006).

Figure 9 shows the histogram of the outflow direction as a function of the angle between the outflow axis and the initial magnetic field direction, for the model with a moderately strong magnetic field (model S1). All the outflows created by the time when SFE ≃ 0.16 are shown. The outflow direction shows a broad distribution, taking its maximum at the mean angle of about 30°. The conclusion is that, although there is some preference for the outflows to align with the global magnetic field direction, the alignment is not as strong as one may naively expect. This is because most of the dense cores are formed through outflow-induced turbulent compression, which can push the magnetic field on the core scale away from the global field direction. On smaller scales, the misalignment in field direction can be further amplified by gravitational collapse. The lack of a complete alignment between the outflow axes and the global field direction is thus expected. It cannot be used to argue against the presence of a moderately strong global magnetic field.

Zoom In Zoom Out Reset image size

Here, we compare the core properties obtained in Section 3 with observations. We use the core data toward two nearby parsec-scale cluster-forming clumps: the ρ Ophiuchus Main Cloud and NGC 1333 in Perseus. In the ρ Ophiuchus region, Maruta et al. (2010) identified 68 H¹³CO⁺ (J = 1–0) cores by applying clumpfind to the three-dimensional position–position–velocity cube data obtained using the Nobeyama 45 m telescope. In NGC 1333, Walsh et al. (2007) identified 93 N₂H⁺ (J = 1–0) dense cores, again using clumpfind. Since both the H¹³CO⁺ (J = 1–0) and N₂H⁺ (J = 1–0) transitions have critical densities of order ∼10⁵ cm⁻³, the H¹³CO⁺ and N₂H⁺ cores should have densities comparable to our simulated cores, which were identified using a threshold density ∼10⁵ cm⁻³. Since the moderately strong field model yields a core velocity dispersion that is in better general agreement with observations than the non-magnetic or weak-field model, we will concentrate on comparing its cores with the ρ Oph and NGC 1333 data.

In Figure 10(a), we compare the one-dimensional FWHM nonthermal velocity-width–radius relations for the ρ Oph cores of Maruta et al. (2010) and our simulated cores (see Section 3). The three-dimensional velocity dispersions discussed in Section 3 were converted into the one-dimensional FWHM velocity widths using dV_NT = (8ln 2/3)^1/2ΔV_NT. Our simulated cores tend to be somewhat smaller than the ρ Oph cores; the average radius of the ρ Oph cores is estimated at 0.04 pc, about 1.7 times that of model S1. Aside from the size difference, the observed and simulated data are almost indistinguishable. Neither shows a clear dependence of linewidth with size, except for a weak trend for the linewidth to increase slightly with radius. Neither follows Larson's linewidth–size relation. The observations are also consistent with the cores formed from initially magnetically subcritical clouds (see Figure 18 of Nakamura & Li 2008). We did not plot the NGC 1333 cores on Figure 10(a) because the core radius in Table 1 of Walsh et al. (2007) was quoted in single significant digit, which artificially introduces a large scatter in the plot. Nevertheless, inspection by eye reveals no clear correlation between the velocity width and core radius, broadly consistent with the data for both the ρ Oph cores and the simulated cores of model S1.

Zoom In Zoom Out Reset image size

Figure 10(b) compares the virial-parameter–mass relations of the simulated and observed cores. The red and blue dots, and black crosses denotes the ρ Oph (Maruta et al. 2010) and NGC 1333 cores (Walsh et al. 2007), and the cores of model S1, respectively. For the NGC 1333 cores, we recalculated the virial masses by including the thermal contribution from the gas with a temperature of 15 K and adopting a distance of 235 pc (Hirota et al. 2008), instead of 300 pc used in Walsh et al. (2007). We also recalculated the LTE masses by adopting the N₂H⁺ fractional abundance relative to H₂ of $X_{\rm N_2H^+} \sim 3\times 10^{-10}$ derived from observations toward ρ Oph because Walsh et al. (2007) determined their adopted N₂H⁺ fractional abundance such that the points are distributed evenly about the line of equality (M_LTE = M_vir) in their M_LTE–M_vir plot. This N₂H⁺ fractional abundance is determined by taking the average of the values obtained toward Oph A (Di Francesco et al. 2004) and B2 subclumps (Friesen et al. 2010). Figure 10(b) indicates that the virial-parameter–core-mass relation of the nearby cluster-forming clumps follows roughly a single power law of α_vir∝M^−2/3_c with a small dispersion, which is again in good agreement with model S1, which has a moderately strong magnetic field. A dynamically important magnetic field may have played a role in core formation in these regions.

Besides the ρ Oph region, Maruta et al. (2010) found that the velocity dispersions of H¹³CO⁺ cores are almost independent of the core radius in Orion A as well, even though the latter is more than a factor of three farther away. The coarser spatial resolution for Orion means that its cores may contain (smaller, ρ Oph-like) cores blending together. If this is true, then the nearly flat velocity–width–radius relation observed in both regions suggests that the inter-core motions among the neighboring cores are almost comparable to the internal motions in the individual cores. Such a feature was pointed out by André et al. (2007) who measured the velocity difference among neighboring cores using N₂H⁺ (J = 1–0). If the velocity–width–radius relation is flat (ΔV∝ constant) and the core mass is proportional to R³_c, then the virial parameter is scaled as α_vir∝R_cΔV²/(GM_c)∝M^−2/3_c, a similar power law to that derived by Bertoldi & McKee (1992). This again suggests the importance of ambient turbulent pressure in dynamics of the cores as discussed in Section 3.

To isolate the effect of magnetic field on clustered star formation, we present in Figure 11(a) the time evolution of SFE for three models with different field strengths but no outflow feedback. Our simulations show that, as expected, the magnetic field tends to slow down star formation, especially at relatively late times. The rate of star formation remains rather high, however, particularly at early times. For example, from the formation of the first star to the time when the SFE reaches ∼16%, the average SFE per global free-fall time star formation rate (SFR_ff; Krumholz & McKee 2005) is estimated to be 29%, 27%, and 21%, respectively, for models N2, W2, and S2. These values are larger than the observationally inferred values of a few percent (Krumholz & Tan 2007). The large SFR_ff in the absence of outflow feedback is in agreement with previous studies (see also the Appendix for the analytic formula of SFR_ff). For example, using MHD smoothed particle hydrodynamics simulations, Price & Bate (2008, 2009) followed the evolution of 50 M_☉ clumps with 0.25 ⩽ β ⩽ ∞ until t ≲ 1.5t_{ff, cl}. Although their initial model clump mass is too small compared to nearby cluster-forming clumps like ρ Oph and NGC 1333 (which have masses of order 10³ M_☉), their star formation efficiencies per free-fall time are estimated to be over 10%. These simulations imply that it is difficult to reduce the SFR to the observed level by a moderately strong magnetic field alone and other factors are needed to significantly retard star formation in cluster-forming clumps.

Zoom In Zoom Out Reset image size

One way to slow down star formation is the feedback from forming stars. Price & Bate (2009) found that radiative feedback from forming stars does not change the global SFE in a parent clump much, although it significantly suppresses the small-scale fragmentation by increasing the temperature in the high-density material near the protostars. Figure 11(b) shows that the SFR is greatly reduced by the inclusion of outflow feedback, even in the presence of a relatively weak magnetic field, in agreement with previous studies (Li & Nakamura 2006; Nakamura & Li 2007; Wang et al. 2010). The average SFE per global free-fall time is estimated to be 8%, 6%, and 4% for models N1 (no magnetic field), W1 (weak field), and S1 (moderately strong field), respectively. For the model with the moderately strong field, the SFR_ff is more or less comparable to the observed value.

The reason for the large reduction in SFR_ff can be seen in Figure 12, where the magnetic energy is plotted against the evolution time, for models W1 and S1. As shown in Figure 12(a), for the moderately strong field case, the total magnetic energy is dominated by the background uniform field that does not contribute to the force balance in the initial cloud. We therefore illustrate in Figure 12(b) the time evolution of the magnetic energy stored in the distorted component that was amplified by supersonic turbulence. Here, we computed the magnetic energy stored in the distorted component by subtracting the initial magnetic energy from the total magnetic energy. For comparison, we also plotted in Figure 12(b) the evolution of the total kinetic energy of the clump. Since the kinetic energy is often dominated by the unbound high-velocity gas associated with the active outflows, we approximate the total kinetic energy of the clump as the kinetic energy of the gas whose velocity is smaller than 10c_s. For both models, the amplified components tend to increase with time. For the model with the moderately strong field, the amplified component begins to oscillate around a level value after a free-fall time. For both models, the magnetic energy of the amplified component becomes comparable to the kinetic energy of the dense gas by the end of the computation, indicating that a quasi-equipartition has been reached (see also Federrath et al. 2011). For the weak field case, the amplified component is more important than the initial uniform field, resulting in a significantly distorted magnetic field structure. Such a highly distorted magnetic field can be seen in the three-dimensional bird's eye view of the density and magnetic field distributions of model W1, which is presented in the left panel of Figure 13. By contrast, the global magnetic field is well ordered for the case of moderately strong initial magnetic field, which is shown in the right panel of Figure 13.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Polarization maps of submillimeter thermal dust emission have recently been obtained for nearby star-forming regions (e.g., Houde et al. 2004; Girart et al. 2009; Matthews et al. 2009). Here, we present the polarization maps derived from the simulation data for the two magnetically supercritical models with different initial magnetic field strengths (models W1 and S1 which have the magnetic flux-to-mass ratios of $\bar{\lambda }=4.3$ and 1.4, respectively). We computed the polarized thermal dust emission from the MHD model following Padoan et al. (2001). We neglect the effect of self-absorption and scattering because we are interested in the thermal dust emission at submillimeter wavelengths. We further assume that the grain properties are constant and the temperature is uniform. The polarization degree is set such that the maximum is equal to 15%. Figure 14 shows the dust polarization maps calculated from models W1 and S1.

Zoom In Zoom Out Reset image size

Only a small portion of the computation box is shown in each panel of Figure 14. As expected from Figure 13, in the presence of a weak magnetic field, the spatial distribution of the polarization vectors has relatively large fluctuations. The polarization degree tends to be smaller in several parts where the magnetic fields are strongly distorted. In the densest parts, the polarization vectors appear more or less parallel to the local elongated structures or dense filaments, whereas there is no clear correlation between the vectors and density distribution in the less dense parts. This may be due to the fact that the local dense filaments are created by turbulent compression that preferentially enhances the magnetic field component transverse to the shock plane. In contrast, in the presence of the strong magnetic field, the filamentary structure is prominent and the filament axes tend to be perpendicular to the polarization vectors that are almost parallel to the initial magnetic field direction, suggesting that the gas flow is channeled preferentially by the strong magnetic field to form the filaments. The polarization observations can thus provide a handle on the magnetic field strength of cluster-forming clumps.

Recent polarization observations of cluster-forming clumps show that the global magnetic field lines are more or less spatially well ordered (e.g., Davis et al. 2000; Houde et al. 2004; Sugitani et al. 2010, 2011). For example, Sugitani et al. (2010) found that the Serpens cloud core is threaded by an hour-glass shaped well-ordered magnetic field and is elongated in the cross-field direction. Sugitani et al. (2011) found that the Serpens South filamentary infrared dark cloud appears to be threaded by more or less straight global magnetic field. The observed spatially ordered magnetic field structures imply that the magnetic fields in the nearby cluster-forming clumps are likely to be at least moderately strong and dynamically important.

In Figures 15(a)–(c), we plot the CMFs for our identified cores for models N1, W1, and S1, respectively. To enlarge the number of sample cores, we added up all the cores identified at the three stages when the SFE has reached 0.08, 0.12, and 0.16. In addition, we added up the cores identified from the runs using a different initial turbulent realization for which different random numbers were used for both the amplitude and phase. The total number of cores so identified is over 500 for each model.

Zoom In Zoom Out Reset image size

The overall shapes of the CMFs are similar in all three cases, suggesting that the shape of the CMF is relatively insensitive to the initial magnetic field strength. This is different from the turbulent fragmentation scenario proposed by Padoan et al. (2001) where a weak magnetic field is required to reproduce the CMF that is similar in shape to the Salpeter IMF (see also Li et al. 2010). In our simulations, a moderately strong magnetic field can produce a CMF that resembles the Salpeter IMF as well.

We note two interesting features of the computed CMFs. First, there appears to be a lack of massive cores compared to the Salpeter stellar IMF. This is not necessarily a problem, because it is unclear whether the nearby cluster-forming clumps such as ρ Oph and NGC 1333 that we aim to simulate would ever produce massive stars. Alternatively, massive stars can grow from initially less massive cores, fed by global collapsing flows toward the bottom of the clump potential (Smith et al. 2008; Wang et al. 2010). Second, the turnover at the low mass end is less clear in the CMFs than in either the Kroupa or Chabrier IMF. This may be because, at the low mass end, most of the cores are gravitationally unbound and may not form stars. Further studies are required to determine the relation between the CMF and IMF.