FORMATION OF MAGNETIZED PRESTELLAR CORES WITH AMBIPOLAR DIFFUSION AND TURBULENCE

CO12 describe a one-dimensional simplified MHD shock system with velocity and magnetic field perpendicular to each other, including a short discussion of oblique shocks. Here we review the oblique shock equations and write them in a more general form to give detailed jump conditions.

We will consider a plane-parallel shock with uniform pre-shock neutral density ρ₀ and ionization-recombination equilibrium everywhere. The shock front is in the x–y plane, the upstream flow is along the z-direction ($\mathbf {v}_{0} = v_{0}\hat{\mathbf {z}}$), and the upstream magnetic field is in the x–z plane, at an angle θ to the inflow ($\mathbf {B}_0 = B_0\sin \theta \hat{\mathbf {x}} + B_0\cos \theta \hat{\mathbf {z}}$) such that B_x = B₀sin θ is the upstream component perpendicular to the flow. The parameters ${\cal M}$ and β (upstream value of the Mach number and plasma parameter) defined in CO12 therefore become

$$\begin{equation} {\cal M} \equiv {\cal M}_z = \frac{v_{0}}{c_{s}},\quad \frac{1}{\beta _0} \equiv \frac{B_{0}^2}{8\pi \rho _{0}c_{s}^2} = \frac{1}{\beta _x} \frac{1}{\sin ^2\theta }. \end{equation} \tag{ 1 }$$

The jump conditions of MHD shocks are described by compression ratios of density and magnetic field:

$$\begin{equation} r_f \equiv \frac{\rho _{n,\ \mathrm{downstream}}}{\rho _{n,\ \mathrm{upstream}}}, \quad r_{B_\perp } \equiv \frac{B_{\perp,\ \mathrm{downstream}}}{B_{\perp,\ \mathrm{upstream}}}. \end{equation} \tag{ 2 }$$

From Equations (A10) and (A14) in CO12, we have

$$\begin{eqnarray} &&\frac{\sin ^2\theta {r_f}^2}{\beta _0}\left(1-\frac{2\cos ^2\theta }{\beta _0{\cal M}^2}\right)^2= \left({\cal M}^2 + 1 + \frac{\sin ^2\theta }{\beta _0} \nonumber \quad- \frac{{\cal M}^2}{r_f} - r_f\right)\left(1-\frac{2 r_f \cos ^2\theta }{\beta _0{\cal M}^2}\right)^2, \end{eqnarray} \tag{ 3 }$$

which can be solved numerically to obtain explicit solution(s) r_{f, exp}(θ). The compression ratio for the magnetic field perpendicular to the inflow is

$$\begin{equation} r_{B_\perp }(\theta) = r_f(\theta)\frac{1 - \frac{2\cos ^2\theta }{\beta _0{\cal M}^2}}{1 - \frac{2 r_f(\theta) \cos ^2\theta }{\beta _0{\cal M}^2}}. \end{equation} \tag{ 4 }$$

Equation (A17) of CO12 gives an analytical approximation to r_f(θ):

$$\begin{equation} r_{f, \mathrm{app}}(\theta) = \frac{\sqrt{\beta _0}{\cal M}}{\sin \theta }\left[\frac{2\sin \theta }{\sqrt{\beta _0}{\cal M}\tan ^2\theta } + \frac{\sqrt{\beta _0}}{2{\cal M}\sin \theta } + 1\right]^{-1}. \end{equation} \tag{ 5 }$$

Since Equation (3) is a quartic function of θ, there are four possible roots of r_f for each angle, and r_f(θ) = const. = 1 (no-shock solution) is always a solution. When θ is large, Equation (3) has one simple root (r_f = 1) and a multiple root with multiplicity =3. When θ drops below a critical value, θ_crit, Equation (3) has four simple roots, which give us four different values of $r_{B_\perp }$. Figure 1 shows the three explicit solutions for r_f and $r_{B_\perp }$ (r_{f, exp}(θ) and r_{B, exp}(θ)), as well as the approximations (r_{f, app}(θ) and r_{B, app}(θ)) that employ Equation (5).

Zoom In Zoom Out Reset image size

The fact that there are multiple solutions for post-shock properties is the consequence of the non-unique Riemann problem in ideal MHD (see discussions in, e.g., Torrilhon 2003; Delmont & Keppens 2011; Takahashi & Yamada 2013), and whether all solutions are physically real is still controversial. The first set of solutions r_{f, exp}(θ), 1 and r_{B, exp}(θ), 1 shown in Figure 1 gives positive r_f and $r_{B_\perp }$, classified as fast MHD shocks (Shu 1992; Draine & McKee 1993), and is the principal oblique shock solution referred to in this contribution.³ The other two solutions for post-shock magnetic field, r_{B, exp}(θ), 2 and r_{B, exp}(θ), 3, both become negative when θ < θ_crit, indicating that the tangential component of the magnetic field to the shock plane is reversed in the post-shock region. These two solutions are commonly specified as intermediate shocks (e.g., Wu 1987; Karimabadi 1995; Inoue & Inutsuka 2007). Among these two field-reversal solutions, we notice that r_{f, exp}(θ), 2 approaches the hydrodynamic jump condition ($r_{f,\mathrm{hydro}} = {\cal M}^2$) when θ → 0, and r_{B, exp}(θ), 2 is smaller in magnitude than other solutions when θ < θ_crit. Thus, we classify this set of solutions r_{f, exp}(θ), 2 and r_{B, exp}(θ), 2 as the quasi-hydrodynamic shock. This quasi-hydrodynamic solution can create gas compression much stronger than the regularly applied fast shock condition and may be the reason that when θ < θ_crit, even ideal MHD simulations can generate shocked layers with relatively high mass-to-flux ratio (see Sections 4 and 5 for more details).

The definition of θ_crit can be derived from Equation (4), which turns negative when $1 - ({2\cos ^2\theta }/{\beta _0{\cal M}^2}) > 0$ and $1 - ({2 r_f(\theta) \cos ^2\theta }/{\beta _0{\cal M}^2}) <0$:

$$\begin{equation} \cos ^2\theta > \cos ^2\theta _\mathrm{crit} = \frac{\beta _0{\cal M}^2}{2 r_f(\theta _\mathrm{crit})}. \end{equation} \tag{ 6 }$$

Using Equation (5) and considering only the terms ${\sim } {\cal M}$, this becomes

$$\begin{equation} \frac{\cos ^2\theta _\mathrm{crit}}{\sin \theta _\mathrm{crit}} \approx \frac{\sqrt{\beta _0} {\cal M}}{2}, \end{equation} \tag{ 7 }$$

or

$$\begin{equation} \sin ^2\theta _\mathrm{crit} + \frac{\sqrt{\beta _0} {\cal M}}{2} \sin \theta _\mathrm{crit} -1 = 0. \end{equation} \tag{ 8 }$$

Assuming θ_crit ≪ 1, this gives

$$\begin{equation} \theta _\mathrm{crit} \sim \frac{2}{\sqrt{\beta _0} {\cal M}} = \sqrt{2}\frac{v_\mathrm{A,0}}{v_0}, \end{equation} \tag{ 9 }$$

where $v_\mathrm{A,0} \equiv B_0/\sqrt{4\pi \rho _0}$ is the Alfvén speed in the cloud. Therefore, the criterion to have multiple solutions, θ < θ_crit, is approximately equivalent to

$$\begin{equation} v_\perp = v_0\sin \theta \lesssim v_0 \cdot \sqrt{2}\frac{v_\mathrm{A,0}}{v_0} \sim v_\mathrm{A,0}\protect, \protect \end{equation} \tag{ 10 }$$

where v_⊥ is the component of the inflow perpendicular to the magnetic field. Though Equation (9) only provides a qualitative approximation⁴ for θ_crit, Equation (10) suggests that when v_⊥/v_{A, 0} is sufficiently small, high-compression quasi-hydrodynamic shocks are possible.

For a core to collapse gravitationally, its self-gravity must overcome both the thermal and magnetic energy. For a given ambient density ρ ≡ μ_nn and assuming spherical symmetry, the mass necessary for gravity to exceed the thermal pressure support (with edge pressure ρc_s²) is the mass of the critical Bonnor-Ebert sphere (see, e.g., Gong & Ostriker 2009):

$$\begin{eqnarray} M_\mathrm{th,sph} &=& 4.18\frac{{c_s}^3}{\sqrt{4\pi G^3 \rho }}\nonumber\\ &=& 4.4 \,M_\odot \left(\frac{T}{10\,\mathrm{K}}\right)^{3/2}\left(\frac{n}{1000\,\mathrm{cm}^{-3}}\right)^{-1/2} \end{eqnarray} \tag{ 11 }$$

(see Section 3.2 for discussion about the value of μ_n). The corresponding length scale at the original ambient density is

$$\begin{eqnarray} R_\mathrm{th,sph} &\equiv & \left(\frac{3 M_\mathrm{th,sph}}{4\pi \rho }\right)^{1/3} = 2.3 \frac{c_s}{\sqrt{4\pi G \rho }}\nonumber\\ &=& 0.26\,\mathrm{pc} \left(\frac{T}{10\,\mathrm{K}}\right)^{1/2}\left(\frac{n}{1000\,\mathrm{cm}^{-3}}\right)^{-1/2}, \end{eqnarray} \tag{ 12 }$$

although the radius of a Bonnor-Ebert sphere with mass given by Equation (11) would be smaller than Equation (12) by 25%, due to internal stratification.

In a magnetized medium with magnetic field B, the ratio of mass to magnetic flux for a region to be magnetically supercritical⁵ can be written as

$$\begin{equation} \frac{M}{\Phi _B}\left |_\mathrm{mag,crit} \equiv \frac{1}{2\pi \sqrt{G}}.\right. \end{equation} \tag{ 13 }$$

With M = 4πR³ρ/3 and Φ_B = πR²B for a spherical volume at ambient density ρ, this gives

$$\begin{eqnarray} M_\mathrm{mag,sph} &=& \frac{9}{128\pi ^2G^{3/2}}\frac{B^3}{\rho ^2}\nonumber\\ &=& 14 \,M_\odot \left(\frac{B}{10 \,\mu \mathrm{G}}\right)^3\left(\frac{n}{1000\,\mathrm{cm}^{-3}}\right)^{-2} \protect \protect \end{eqnarray} \tag{ 14 }$$

and

$$\begin{equation} R_\mathrm{mag,sph} = \frac{3}{8\pi \sqrt{G}}\frac{B}{\rho } = 0.4\,\mathrm{pc} \left(\frac{B}{10 \,\mu \mathrm{G}}\right) \left(\frac{n}{1000\,\mathrm{cm}^{-3}}\right)^{-1}\protect. \protect \protect \protect \end{equation} \tag{ 15 }$$

A spherical region must have M > M_{th, sph} as well as M > M_{mag, sph} to be able to collapse. In the cloud environment (the pre-shock region), B ∼ 10 μG and n ∼ 1000 cm⁻³ are typical. Comparing Equations (11) and (14), the magnetic condition is more strict than the thermal condition; if cores formed from a spherical volume, the mass would have to exceed ∼10 M_☉ in order to collapse. This value is much larger than the typical core mass (∼1 M_☉) identified in observations. This discrepancy is the reason why traditionally ambipolar diffusion is invoked to explain how low-mass cores become supercritical.

We can examine the ability for magnetically supercritical cores to form isotropically in a post-shock layer. The normalized mass-to-flux ratio

$$\begin{equation} \Gamma \equiv \frac{M}{\Phi _B}\cdot 2\pi \sqrt{G} \end{equation} \tag{ 16 }$$

of a spherical volume with density ρ, magnetic field B, and mass M is

$$\begin{eqnarray} \Gamma _\mathrm{sph} &=& \frac{8\pi \sqrt{G}}{3}\left(\frac{3}{4\pi }\right)^{1/3} M^{1/3}\rho ^{2/3}B^{-1}\nonumber \\ &=& 0.4 \left(\frac{M}{M_\odot }\right)^{1/3}\left(\frac{n}{1000\,\mathrm{cm}^{-3}}\right)^{2/3}\left(\frac{B}{10 \,\mu \mathrm{G}}\right)^{-1}.\quad\quad \end{eqnarray} \tag{ 17 }$$

Or, with Σ = 4Rρ/3 ≡ μ_nN_n for a sphere, we have

$$\begin{equation} \Gamma _\mathrm{sph} =2\pi \sqrt{G}\cdot \frac{\Sigma }{B} = 0.6\left(\frac{N_n}{10^{21}\,\mathrm{cm}^{-2}}\right)\left(\frac{B}{10 \,\mu \mathrm{G}}\right)^{-1}. \end{equation} \tag{ 18 }$$

Considering the cloud parameters from Figure 1 (${\cal M} = 10$, B₀ = 10 μG, n₀ = 1000 cm⁻³), the post-shock density and magnetic field are approximately n_ps ∼ 10⁴ cm⁻³ and B_ps ∼ 50 μG when θ > θ_crit. A solar-mass spherical region in this shocked layer will have Γ_{ps, sph} ≈ 0.37; spherical contraction induced by gravity would be suppressed by magnetic fields. Thus, typical post-shock conditions are unfavorable for forming low-mass cores by spherical contraction in ideal MHD.

Furthermore, using r_f and $r_{B_\perp }$ defined in Section 2.1, we can compare Γ_{ps, sph} and the pre-shock value Γ_{pre, sph} for spherical post-shock and pre-shock regions:

$$\begin{eqnarray} \frac{\Gamma _\mathrm{ps,sph}}{\Gamma _\mathrm{pre,sph}} &=& \left(\frac{M_\mathrm{ps}}{M_\mathrm{pre}}\right)^{1/3}\left(\frac{\rho _\mathrm{ps}}{\rho _\mathrm{pre}}\right)^{2/3}\left(\frac{B_\mathrm{ps}}{B_\mathrm{pre}}\right)^{-1}\nonumber\\ &\approx & \left(\frac{M_\mathrm{ps}}{M_\mathrm{pre}}\right)^{1/3} {r_f}^{2/3} {r_{B_\perp }}^{-1}. \end{eqnarray} \tag{ 19 }$$

Considering volumes containing similar mass, M_ps ∼ M_pre, the ratio between the post-shock and pre-shock Γ_sph is smaller than unity when θ > θ_crit, because Equation (4) shows that $r_{B_\perp }$ is larger than r_f. Thus, provided that θ > θ_crit, the post-shock layer will actually have stronger magnetic support than the pre-shock region for a given spherical mass.

Based on the above considerations, formation of low-mass supercritical cores appears difficult in ideal MHD. Adapting classical ideas, one might imagine that low-mass subcritical cores form quasi-statically within the post-shock layer and then gradually lose magnetic support via ambipolar diffusion to become magnetically supercritical in a timescale ∼1–10 Myr. A process of this kind would, however, give prestellar core lifetimes longer than observed, and most cores would have Γ < 1 (inconsistent with observations).

Two alternative scenarios could lead to supercritical core formation in a turbulent magnetized medium. First, the dynamic effects during a turbulence-induced shock (including rapid, transient ambipolar diffusion and the quasi-hydrodynamic compression when θ < θ_crit) may increase the compression ratio of neutrals, creating $r_f \gg r_{B_\perp }$ and Γ_{ps, sph} > 1, enabling low-mass supercritical cores to form. Second, even if the post-shock region is strongly magnetized, mass can accumulate through anisotropic condensation along the magnetic field until both the thermal and magnetic criteria are simultaneously satisfied. In this study, we carefully investigate these two scenarios, showing that both effects contribute to the formation of low-mass supercritical cores within timescale ≲ 0.6 Myr, regardless of ionization or magnetic obliquity.

To examine core formation in shocked layers of partially ionized gas, we employ a three-dimensional convergent flow model with ambipolar diffusion, self-gravity, and a perturbed turbulent velocity field. We conducted our numerical simulations using the Athena MHD code (Stone et al. 2008) with Roe's Riemann solver. To avoid negative densities if the second-order solution fails, we instead use first-order fluxes for bad zones. The self-gravity of the domain, with an open boundary in one direction and periodic boundaries in the other two, is calculated using the fast Fourier transformation method developed by Koyama & Ostriker (2009). Ambipolar diffusion is treated in the strong coupling approximation, as described in Bai & Stone (2011), with super time-stepping (Choi et al. 2009) to accelerate the evolution.

The equations we solve are

$$\begin{eqnarray} \frac{\partial \rho _n}{\partial t} &+ \mathbf {\nabla }\cdot \left(\rho _n\mathbf {v}\right) = 0, \end{eqnarray} \tag{ 20a }$$

$$\begin{eqnarray} \frac{\partial \rho _n\mathbf {v}}{\partial t} &+ \mathbf {\nabla }\cdot \left(\rho _n\mathbf {v}\mathbf {v} - \frac{\mathbf {B}\mathbf {B}}{4\pi }\right) + \mathbf {\nabla }P^* = 0, \end{eqnarray} \tag{ 20b }$$

$$\begin{eqnarray} \frac{\partial \mathbf {B}}{\partial t} &+ \mathbf {\nabla }\times \left(\mathbf {B}\times \mathbf {v}\right) = \mathbf {\nabla }\times \left[\frac{\left(\left(\mathbf {\nabla }\times \mathbf {B}\right)\times \mathbf {B}\right)\times \mathbf {B}}{4\pi \rho _i\rho _n\alpha }\right],\quad\quad \end{eqnarray} \tag{ 20c }$$

where P* = P + B²/(8π). For simplicity, we adopt an isothermal equation of state P = ρc_s². The numerical setup for inflow and turbulence is similar to that adopted by GO11. For both the whole simulation box initially and the inflowing gas subsequently, we apply perturbations following a Gaussian random distribution with a Fourier power spectrum as described in GO11. The scaling law for supersonic turbulence in GMCs obeys the relation

$$\begin{equation} \frac{\delta v_\mathrm{1D} (\ell)}{\sigma _{v,\mathrm{cloud}}} = \left(\frac{\ell }{2R_\mathrm{cloud}}\right)^{1/2}, \end{equation} \tag{ 21 }$$

where δv_1D(ℓ) represents the one-dimensional velocity dispersion at scale ℓ, and σ_{v, cloud} is the cloud-scale one-dimensional velocity dispersion. In terms of the virial parameter α_vir ≡ 5σ_v²R_cloud/(GM_cloud) with M_cloud ≡ 4πρ₀R_cloud³/3, and for the inflow Mach number ${\cal M}$ comparable to σ_v/c_s of the whole cloud, the three-dimensional velocity dispersion $\delta v=\sqrt{3}\cdot \delta v_\mathrm{1D}$ at the scale of the simulation box would be

$$\begin{equation} \delta v (L_\mathrm{box}) = \sqrt{3} \left(\frac{\pi G \alpha _\mathrm{vir}}{15}\right)^{1/4} {\cal M}^{1/2} {c_s}^{1/2} {\rho _0}^{1/4} {L_\mathrm{box}}^{1/2}. \end{equation} \tag{ 22 }$$

To emphasize the influence of the cloud magnetization instead of the perturbation field, our simulations are conducted with 10% of the value δv(L_box), or δv = 0.14 km s⁻¹ with α_vir = 2. With larger δv(L_box), simulations can still form cores, but because non-self-gravitating clumps can easily be destroyed by strong velocity perturbations and no core can form before the turbulent energy dissipates, it takes much longer, with corresponding higher computational expense.

A schematic showing our model setup is shown in Figure 2. Our simulation box is 1 pc on each side and represents a region within a GMC where a large-scale supersonic converging flow with velocity v₀ and −v₀ (i.e., in the center-of-momentum frame) collides. The z-direction is the large-scale inflow direction, and we adopt periodic boundary conditions in the x- and y-directions. We initialize the background magnetic field in the cloud, B₀, in the x–z plane, with an angle θ with respect to the convergent flow. For simplicity, we treat the gas as isothermal at temperature T = 10 K, such that the sound speed is c_s = 0.2 km s⁻¹. The neutral density within the cloud, ρ₀, is set to be uniform in the initial conditions and in the upstream converging flow.

Zoom In Zoom Out Reset image size

It has been shown that ionization-recombination equilibrium generally provides a good approximation to the ionization fraction within GMCs for the regime under investigation (CO12). Thus, the number density of ions in our model can be written as

$$\begin{equation} n_i = \frac{\rho _i}{\mu _i}= 10^{-6} \chi _{i0} \left(\frac{\rho _n}{\mu _n}\right)^{1/2}, \end{equation} \tag{ 23 }$$

with

$$\begin{equation} \chi _{i0} \equiv 10^6 \times \sqrt{\frac{\zeta _\mathrm{CR}}{\alpha _\mathrm{gas}}} \end{equation} \tag{ 24 }$$

determined by the cosmic-ray ionization rate (ζ_CR) and the gas-phase recombination rate (α_gas). The ionization coefficient, χ_i0, has values ∼1–20 (McKee et al. 2010) and is the model parameter that controls ambipolar diffusion effects in our simulations, following CO12. We use typical values of the mean neutral and ion molecular weight μ_n and μ_i of 2.3m_H and 30m_H, respectively, which give the collision coefficient (see Equation (20c)) between neutrals and ions α = 3.7 × 10¹³ cm³ s⁻¹ g⁻¹.

The physical parameters defining each model are ρ₀, v₀ = |v₀|, B₀ = |B₀|, θ, and χ_i0. We set the upstream neutral number density to be n₀ = ρ₀/μ_n = 1000 cm⁻³ in all simulations, consistent with typical mean molecular densities within GMCs⁶ (e.g., Larson 1981; Williams et al. 2000; Bot et al. 2007; Bolatto et al. 2008). We choose the upstream B₀ = 10 μG as typical of GMC values (Goodman et al. 1989; Crutcher et al. 1993; Heiles & Crutcher 2005; Heiles & Troland 2005) for all our simulations. To keep the total number of simulations practical, we set the large-scale inflow Mach number to ${\cal M}=10$ for all models. Exploration of the dependence on Mach number of ambipolar diffusion and of core formation has been studied in previous simulations (CO12 and GO11, respectively). For our parameter survey, we choose θ = 5, 20, and 45 degrees to represent small (θ < θ_crit), intermediate (θ > θ_crit), and large (θ ≫ θ_crit) angles between the inflow velocity and cloud magnetic field. For each θ, we conduct simulations with χ_i0 = 3, 10, and ideal MHD to cover situations with strong, weak, and no ambipolar diffusion. We also run corresponding hydrodynamic simulations with the same ρ₀ and v₀ for comparison.

A full list of models is contained in Table 1. Table 1 also lists the steady-state post-shock properties, as described in Section 2.1. Solutions for all three types of shocks are listed for the θ = 5° (A5) case. For the θ = 20° and θ = 45° cases, there is only one shock solution. Also included in Table 1 are the nominal values of critical mass and radius for spherically symmetric volumes to be self-gravitating under these steady-state post-shock conditions, as discussed in Section 2.2 (see Equations (11), (12), (14), and (15)). Both "thermal" and "magnetic" critical masses are listed. In most models, M_{mag, sph} > M_{th, sph} and M_{mag, sph} ≫ M_☉, indicating that the post-shock regions are dominated by magnetic support, and either ambipolar diffusion or anisotropic condensation would be needed to form low-mass supercritical cores, as discussed in Section 2.2. On the other hand, the quasi-hydrodynamic shock solution for models with θ < θ_crit (i.e., A5 cases) has M_{mag, sph} < M_{th, sph} < M_☉ downstream. If this shock solution could be sustained, then in principle low-mass supercritical cores could form by spherical condensation of post-shock gas.

In order to collect sufficient statistical information on the core properties from simulations, we repeat each parameter set six times with different random realizations of the same perturbation power spectrum for the turbulence. The resolution is 256³ for all simulations such that Δx ≈ 0.004 pc, or ∼800 AU. We tested this setup with two times of this resolution (Δx ≈ 0.002 pc), and the resulting dense structures are highly similar. Though the individual core properties vary around ±50%, the median values (which are more important in our statistical study) only change within ±10%–30%. Thus, our simulations with Δx ≈ 0.004 pc are well resolved for investigations of core properties.

To measure the physical properties of the cores formed in our simulations, we apply the GRID core-finding method developed by GO11, which uses gravitational potential isosurfaces to identify cores. In this approach, the largest closed potential contour around a single local minimum of the gravitational potential defines the material eligible to be part of a core. We define the bound core region as all the material within the largest closed contour that has the sum of gravitational, magnetic, and thermal energy negative.⁷ All of our cores are, by definition, self-gravitating.

The essential quantity to measure the significance of magnetic fields in self-gravitating cores is the ratio of mass to magnetic flux (Mestel & Spitzer 1956; Mouschovias & Spitzer 1976). From Gauss's law the net flux of the magnetic field through a closed surface is always zero. As a result, to measure the magnetic flux within a core, we need to first define a cross section of the core and then measure the net magnetic flux through the surface of the core defined by this cross section (which is the same as the flux through the cross section itself).

To define the cross section through a core, we use the plane perpendicular to the average magnetic field that also includes the minimum of the core's gravitational potential. This choice ensures that we measure the magnetic flux through the part of the core with strongest gravity. After defining this plane, we separate the core into an upper half and a lower half and measure the magnetic flux Φ_B through one of the halves. In practice, we compute this by first finding all zones that contain at least one face that is on the core surface, and we assign normal vectors $\mathbf {{\hat{n}}}$ (pointing outward) to those faces. From these, we select only those in the upper "hemisphere" of the core. After we have a complete set of those grid faces that are on the upper half of the core surface, we sum up their $\mathbf {B}\cdot \mathbf {{\hat{n}}}$ to get the net magnetic flux of the core. This method is tested in spherical and rectangular "cores" with magnetic fields in arbitrary directions. Note that this method works best when the core is approximately spherical (without corners).

After we have the measurement of magnetic flux Φ_B, we can calculate the mass-to-flux ratio of the core, M/Φ_B. This determines whether the magnetic field can support a cloud against its own self-gravity. The critical value of M/Φ_B differs with the geometry of the cloud, but the value varies only within ∼10% (e.g., Mouschovias & Spitzer 1976; Nakano & Nakamura 1978; Tomisaka et al. 1988, or see review in McKee & Ostriker 2007). We therefore choose the commonly used value $(2\pi \sqrt{G})^{-1}$ (e.g., Kudoh & Basu 2011; Vázquez-Semadeni et al. 2011; CO12) as a reference value and define the normalized mass-to-flux ratio as $\Gamma \equiv 2\pi \sqrt{G}\cdot M/\Phi _B$ (see Equation (16)). For a prestellar core with Γ > 1, the gravitational force exceeds the magnetic support and the core is magnetically supercritical. A subcritical core has Γ < 1 and is ineligible for wholesale collapse unless magnetic fields diffuse out.

Figure 8 shows the distribution of mass and size of cores for all model parameters. The masses range from 0.04 to 2.5 M_☉, with a peak around ∼0.6 M_☉; the core sizes are between 0.015 and 0.07 pc, with a peak around ∼0.03 pc. These are consistent with observational results (e.g., Motte et al. 2001; Ikeda et al. 2009; Rathborne et al. 2009; Kirk et al. 2013). Also, the distribution of the core mass shows a similar shape to the observed core mass function (e.g., Simpson et al. 2008; Rathborne et al. 2009; Könyves et al. 2010). Interestingly, the peak in the distribution is close to the value given by Equation (7) from GO11:

$$\begin{equation} M_\mathrm{BE,\ ps} = 1.2 \frac{{c_s}^4}{\sqrt{G^3 P_\mathrm{ps}}}= 1.2 \frac{{c_s}^3}{\sqrt{G^3 \rho _0}}\frac{1}{{\cal M}} \rightarrow 0.45\,M_\odot. \end{equation} \tag{ 27 }$$

This mass is characteristic of what is expected for collapse of a thermally supported core that is confined by an ambient medium with pressure equal to the post-shock value,⁹ where the numerical figure uses values for the mean cloud density and large-scale Mach number equal to those of the converging flow in our simulations, n₀ = 1000 cm⁻³ and ${\cal M} = 10$. Correspondingly, since the critical ratio of mass and radius is M_BE/R_BE = 2.4c_s²/G (Bonnor 1956), the characteristic size expected for a collapsing core formed in a post-shock region when the Mach number of the large-scale converging flow is ${\cal M}$ and the mean cloud density is ρ₀ is

$$\begin{equation} L_\mathrm{BE} = 2 R_\mathrm{BE} = \frac{c_s}{\sqrt{G\rho _0}}\frac{1}{{\cal M}} \rightarrow 0.04\,\mathrm{pc}. \end{equation} \tag{ 28 }$$

This is again comparable to the peak value of the core size distribution in Figure 8.

Zoom In Zoom Out Reset image size

We also separately explore the dependence of core mass, size, magnetic field strength, and mass-to-flux ratio on model parameters, as shown in Figure 9. Our results show that the core mass is relatively insensitive to both the ionization (i.e., ambipolar diffusion effect) and obliquity of the upstream magnetic field (Figure 9, top left). The median masses are within a factor 2.4 of the mean of the whole distribution, 0.68 M_☉, or a factor two of the median of all core masses (0.47 M_☉). Similarly, median core sizes vary only between 0.022 pc and 0.034 pc for the various parameter sets, with a median of 0.03 pc. Note that we chose to compare median values between different parameter sets in Figure 9 instead of mean values used in Table 3, because an average can be affected by any single value being high or low relative to the other samples. The median value, on the other hand, represents the central tendency better, and with the ±25% values we can have a better understanding of the sample distribution.

Zoom In Zoom Out Reset image size

We note in particular that for the θ = 20° and θ = 45° ideal MHD cases, the masses in Figures 8 and 9 are more than an order of magnitude lower than the limits for a spherical region at post-shock conditions to be magnetically supercritical, as listed in Tables 1 and 2. This implies that the low-mass bound cores found in the simulations did not form isotropically. We discuss this further in Section 6.

To further investigate the relationship between core masses and sizes, we binned the data set by log L_core and calculate the average core mass and mean density for different model parameters. The results are shown in Figure 10, where we chose four models with different magnetization and ionization levels to compare: HD (hydrodynamics; no magnetization), A5X3 (low ionization, weak upstream magnetic field parallel to the shock), A20X10 (moderate ionization and magnetic field), and A45ID (ideal MHD, strong magnetic field). In both the mass-size and density-size plots, the differences among models are small, and all four curves have similar shape. In fact, from all resolved cores identified in our simulations, we found a power-law relationship between the core mass and size, M∝L^k, with best-fitted value k = 2.28. This is consistent with many core-property surveys toward different molecular clouds (e.g., Elmegreen & Falgarone 1996; Curtis & Richer 2010; Roman-Duval et al. 2010; Kirk et al. 2013), in which k = 1.2–2.4 with various molecule tracers (for more details, see Figure 7 and corresponding discussions in Kirk et al. 2013).

Zoom In Zoom Out Reset image size

Figure 11 shows the distribution of core mass-to-flux ratio, a roughly normal distribution with peak at Γ ∼ 3. This range of Γ is quite similar to observational results (Γ ∼ 1–4; Falgarone et al. 2008; Troland & Crutcher 2008). In addition, the color-coded histogram in Figure 11 shows how the mass-to-flux ratio depends on magnetization: the high-end region (Γ ≳ 5) is composed of blue-green pieces (which represent models with lower ionization), while the low-end tail is mostly red and orange (highly ionized models). Note that essentially all of the cores in our simulations are magnetically supercritical (Γ > 1), which is self-consistent with our core-finding criterion (gravitationally bound; E_g + E_th + E_B < 0).

Zoom In Zoom Out Reset image size

The tendency of models with lower ionization to form cores with higher mass-to-flux ratio is very clearly seen in Figure 9 (bottom right). The median value of the core mass-to-flux ratio also decreases with increasing θ as the value of the upstream B_x = B₀sin θ increases. Also from Figure 9 (top right), the average core magnetic fields show a similar tendency as in post-shock magnetic field (see Table 2), which decrease at lower ionization fractions for models with the same pre-shock magnetic field structure (same θ). The larger and more systematic variation of $\overline{B}$ than M with model parameters suggests that the core mass-to-flux ratio is not decided by the core mass, but by the core magnetic field. This is also shown in Figure 12, where we binned the data by M_core and calculated the average core mass-to-flux ratio in each bin for different models. For cores with similar mass, the mass-to-flux ratios of cores formed in environments with low ionization and magnetization are much higher than those with stronger and better coupled magnetic fields.

Zoom In Zoom Out Reset image size

The fact that the median value of magnetic field strength within the core depends on pre-shock magnetic obliquity and ionization is consistent with our discussions in Section 4 that magnetic fields are lower in shocked regions that have longer transient timescales. Since lower ionization fraction leads to stronger ambipolar diffusion and a longer transient stage,¹⁰ it is logical to expect that the cores formed in weakly ionized clouds have lower magnetic field than those formed with higher ionization fraction (or strongly coupled ions and neutrals).

In addition, Figure 9 (top right) shows that cores formed in models with small θ (A5 cases) have weaker magnetic fields inside even with higher ionization fraction or ideal MHD, which indicates that the magnetic field is less compressed by the shock when the inflow is almost parallel to the upstream magnetic field. This is consistent with the discussion in Section 2: when θ < θ_crit, the MHD shock becomes a composite compounded of the regular (fast) mode and the quasi-hydrodynamic mode, which has a relatively small magnetic field compression ratio. Thus, the magnetic obliquity relative to the shock has a similar effect to the cloud ionization fraction in determining field strengths in prestellar cores.

Based on the results shown in Sections 5.1 and 5.2, we conclude that magnetic effects do not appear to control core mass and size. This suggests that once a core becomes strongly gravitationally bound, magnetic effects are relatively unimportant to its internal structure. However, the formation process of gravitationally bound cores is highly dependent on magnetic effects. As noted above, Figure 6 shows clear differences in the large-scale structures from which cores condense; we discuss core formation further in Section 6. Also, cores are born with either lower or higher magnetic field, depending on the magnetic field structure and the ambipolar diffusion in their surrounding environment.

The fact that gravitationally supercritical low-mass cores (with M ≪ M_{mag, sph}) can form in the highly magnetized post-shock medium even without ambipolar diffusion suggests that these cores did not contract isotropically. Figure 13 provides a close-up view of the core forming process in a highly magnetized environment with ideal MHD, from model A20ID. At stages earlier than shown, the directions of the perturbed magnetic field and gas velocity are determined randomly by the local turbulence. The magnetic field is compressed by the shock (similar to Figures 3−5), such that in the post-shock layer it is nearly parallel to the shock front (along $\hat{\mathbf {x}}$). When the magnetic field strength increases, the velocity is forced to become increasingly aligned parallel to the flow, as shown in Figure 13. By the time t = 0.65 Myr, a very dense core has formed by gathering material preferentially along the magnetic field lines. After the core becomes sufficiently massive, its self-gravity will distort the magnetic field and drag material inward even in the direction perpendicular to the magnetic field lines (t = 0.77 Myr; Figure 13). This collapsing process with a preferential direction is similar to the post-shock focusing flows found in previous studies (e.g., Inoue & Fukui 2013; Vaidya et al. 2013) where the gas is confined by the strong magnetic field in the shock-compressed region.

Zoom In Zoom Out Reset image size

In fact, anisotropic condensation is key to core formation not only with ideal MHD, but for all cases. Figure 14 shows space-time diagrams of all three velocity components (v_x, v_y, v_z) around collapsing cores in different parameter sets. Though models with stronger post-shock magnetic fields (larger θ, larger χ_i0, or ideal MHD) have more dominant v_x (bluer/redder in the color map), there is a general preference to condense preferentially along the magnetic field lines (in the x-direction) among all models, regardless of upstream magnetic obliquity and the ambipolar diffusion level. Figure 14 also shows that there are flows perpendicular to the mean magnetic field (along $\mathbf {\hat{y}}$ or $\mathbf {\hat{z}}$) in the final ∼0.1 Myr of the simulations, indicating the stage of core collapse. The prominent gas movement along $\mathbf {\hat{x}}$ long before each core starts to collapse shows that cores acquire masses anisotropically along the magnetic field lines, and thus anisotropic condensation is important for all models.

Zoom In Zoom Out Reset image size

We have shown in Section 2.2 quantitatively that isotropic formation of low-mass supercritical cores is not possible for oblique shocks with ideal MHD, because the minimum mass for a spherical volume to become magnetically supercritical is ⩾10 M_☉ (see Equation (14) and the M_{mag, sph} entries in Tables 1 and 2 for cases A20ID and A45ID), much larger than the typical core mass (∼1 M_☉). However, non-spherical regions may have smaller critical mass. Consider, for example, a core that originates as a prolate spheroid with semi-major axis R₁ along the magnetic field and semi-minor axis R₂ perpendicular. The mass-to-flux ratio is then

$$\begin{equation} \frac{M}{\Phi _B}\left |_\mathrm{prolate} = \frac{4\pi R_1 {R_2}^2 \rho /3}{\pi {R_2}^2 B} = \frac{4}{3}\frac{R_1 \rho }{B}.\right. \end{equation} \tag{ 29 }$$

The critical value for R₁ would be the same as given in Equation (15), but the critical mass would be lower than that in Equation (14) by a factor (R₂/R₁)². For R₁/R₂ ∼ 3–4, the critical mass will be similar to that found in the simulations.

Here we provide a physical picture for core formation via initial flow along the magnetic field, as illustrated in Figures 13 and 14. Consider a post-shock layer with density ρ_ps and magnetic field B_ps. For a cylinder with length L along the magnetic field and radius R, the normalized mass-to-flux ratio is

$$\begin{equation} \Gamma _\mathrm{cyl} = \frac{\pi R^2 L\rho _\mathrm{ps}}{\pi R^2 B_\mathrm{ps}}\cdot 2\pi \sqrt{G}, \end{equation} \tag{ 30 }$$

and the critical length along the magnetic field for it to be supercritical is

$$\begin{equation} L_\mathrm{mag, cyl} = \frac{B_\mathrm{ps}}{\rho _\mathrm{ps}}\frac{1}{2\pi \sqrt{G}} \end{equation} \tag{ 31 }$$

(note that up to a factor 3/4, this is the same as Equation (15)). The critical mass M_{mag, cyl} = πR²L_{mag, cyl}ρ_ps can then be written as

$$\begin{equation} M_\mathrm{mag,cyl} = \frac{R^2B_\mathrm{ps}}{2\sqrt{G}} = 1.2 \,M_\odot \left(\frac{R}{0.05\,\mathrm{pc}}\right)^2\left(\frac{B_\mathrm{ps}}{50 \,\mu \mathrm{G}}\right). \end{equation} \tag{ 32 }$$

This cylinder is gravitationally stable to transverse contraction unless L ≲ 2R (Mestel & Spitzer 1956). However, contraction along the length of the cylinder is unimpeded by the magnetic field and will be able to overcome pressure forces provided that L_{mag, cyl} exceeds the thermal Jeans length, which is true in general for oblique shocks in typical conditions under consideration. The longitudinal contraction will produce an approximately isotropic core of radius R when the density has increased by a factor

$$\begin{equation} \frac{\rho ^{\prime }}{\rho _\mathrm{ps}} = \frac{L_\mathrm{mag,cyl}}{2R}, \end{equation} \tag{ 33 }$$

and at this point transverse contraction would no longer be magnetically impeded. For the core to have sufficient self-gravity to overcome thermal pressure at this point, the radius would have to be comparable to R_{th, sph} (see Equation (12)):

$$\begin{equation} R \sim R_\mathrm{th,sph} = 2.3\frac{c_s}{\sqrt{4\pi G \rho ^{\prime }}}. \end{equation} \tag{ 34 }$$

Combining Equations (31), (33), and (34) yields

$$\begin{equation} \rho ^{\prime } = 0.19 \frac{{B_\mathrm{ps}}^2}{4\pi {c_s}^2} \protect \protect \end{equation} \tag{ 35 }$$

and

$$\begin{equation} R = 5.3\frac{{c_s}^2}{\sqrt{G}B_\mathrm{ps}} = 0.05\,\mathrm{pc} \left(\frac{B_\mathrm{ps}}{50 \,\mu \mathrm{G}}\right)^{-1}\left(\frac{T}{10\,\mathrm{K}}\right). \end{equation} \tag{ 36 }$$

Substituting Equation (36) in Equation (32), the minimum mass that will be both magnetically and thermally supercritical, allowing for anisotropic condensation along B, will be

$$\begin{equation} M_\mathrm{crit} = 14 \frac{{c_s}^4}{G^{3/2}B_\mathrm{ps}} = 1.3\,M_\odot \left(\frac{B_\mathrm{ps}}{50\,\mu \mathrm{G}}\right)^{-1}\left(\frac{T}{10\,\mathrm{K}}\right)^2. \end{equation} \tag{ 37 }$$

Thus, anisotropic contraction can lead to low-mass supercritical cores, with values comparable to those formed in our simulations.¹¹

In addition, anisotropic condensation also helps explain why the core masses are quite similar for HD and MHD models, and independent of the angle between upstream magnetic field and converging flow. Note that Equation (37) only depends on the post-shock magnetic field strength. For a magnetized shock, the post-shock magnetic pressure must balance the pre-shock momentum flux: B_ps²/8π ∼ ρ₀v₀². Therefore, Equation (37) can be expressed as

$$\begin{eqnarray} M_\mathrm{crit} &=& 2.8 \frac{{c_s}^4}{\sqrt{G^3 \rho _0 {v_0}^2}} = 2.1\,M_\odot \left(\frac{n_0}{1000\,\mathrm{cm}^{-3}}\right)^{-1/2}\nonumber\\ &&\times\left(\frac{v_0}{1\,\mathrm{km\,s^{-1}}}\right)^{-1} \left(\frac{T}{10\,\mathrm{K}}\right)^2. \end{eqnarray} \tag{ 38 }$$

This is equivalent to Equation (24) of GO11 with ψ = 2.8. GO11 also pointed out that ρ₀v₀² will be proportional to GΣ_GMC² for a gravitationally bound turbulence-supported GMC. Thus, using Equation (28) of GO11 in a cloud with virial parameter α_vir, Equation (38) would become

$$\begin{equation} M_\mathrm{crit}=2.8 \,M_\odot \left(\frac{T}{10\,\mathrm{K}}\right)^2\left(\frac{\Sigma _\mathrm{GMC}}{100 \,M_\odot \,\mathrm{pc}^{-2}}\right)^{-1}{\alpha _\mathrm{vir}}^{-1/2}. \end{equation} \tag{ 39 }$$

Equations (38) and (39) suggest that M_crit is not just independent of magnetic field direction upstream, but also independent of magnetic field strength upstream. That is, when cores form in post-shock regions (assuming that the GMC is magnetically supercritical at large scales), the critical mass is determined by the dynamical pressure in the cloud, independent of the cloud's magnetization. The models studied here all have the same dynamical pressure ρ₀v₀² and same upstream B₀. It will be very interesting to test whether for varying B₀ the core masses remain the same, and whether the scaling proposed in Equation (38) holds for varying total dynamic pressure.