FORMATION OF COLLAPSING CORES IN SUBCRITICAL MAGNETIC CLOUDS: THREE-DIMENSIONAL MAGNETOHYDRODYNAMIC SIMULATIONS WITH AMBIPOLAR DIFFUSION

We solve the three-dimensional MHD equations including self-gravity and ambipolar diffusion, assuming that neutrals are much more numerous than ions:

$$\begin{equation} \frac{\partial \rho }{\partial t} + \mbox{\boldmath $v$} \cdot \nabla \rho = -\rho \nabla \cdot \mbox{\boldmath $v$}, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial \mbox{\boldmath $v$}}{\partial t} + (\mbox{\boldmath $v$} \cdot \nabla) \mbox{\boldmath $v$} =-\frac{1}{\rho } \nabla p + \frac{1}{c \rho } \mbox{\boldmath $j$} \times \mbox{\boldmath $B$} - \nabla \psi, \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{\partial \mbox{\boldmath $B$}}{\partial t} = \nabla \times (\mbox{\boldmath $v$} \times \mbox{\boldmath $B$}) + \nabla \times \left[\frac{\tau _{ni}}{c\rho } (\mbox{\boldmath $j$} \times \mbox{\boldmath $B$}) \times \mbox{\boldmath $B$}\right], \end{equation} \tag{ 3 }$$

$$\begin{equation} \mbox{\boldmath $j$}=\frac{c}{4\pi } \nabla \times \mbox{\boldmath $B$}, \end{equation} \tag{ 4 }$$

$$\begin{equation} \nabla ^2 \psi = 4\pi G \rho, \end{equation} \tag{ 5 }$$

$$\begin{equation} p=c_s^2 \rho, \end{equation} \tag{ 6 }$$

where ρ is the density of neutral gas, p is the pressure, v is the velocity, B is the magnetic field, j is the electric current density, ψ is the self-gravitating potential, and c_s is the sound speed. Here, we consider the molecular clouds whose number density is about 10³–10⁶ cm⁻³. In this case, the cooling time is generally much shorter than the dynamical time. Therefore, instead of solving a detailed energy equation, we assume isothermality for each Lagrangian fluid particle (Kudoh & Basu 2003, 2006):

$$\begin{equation} \frac{d c_s}{dt} = \frac{\partial c_s}{\partial t} + \mbox{\boldmath $v$} \cdot \nabla c_s = 0. \end{equation} \tag{ 7 }$$

This means that each parcel of molecular cloud gas as well as surrounding warm gas (see Section 2.2) retains its initial temperature. For the neutral–ion collision time in Equation (3) and associated quantities, we follow Basu & Mouschovias (1994), so that

$$\begin{equation} \tau _{\rm ni} = 1.4\, \frac{m_{\rm i}+m_{\rm n}}{\rho _{\rm i} \langle \sigma w \rangle _{\rm in}}, \end{equation} \tag{ 8 }$$

where ρ_i is the density of ions and 〈σw〉_in is the average collisional rate between ions of mass m_i and neutrals of mass m_n. Here, we use typical values of HCO⁺–H₂ collisions, for which 〈σw〉_in = 1.69 × 10⁻⁹ cm⁻³ s⁻¹, m_i/m_n = 14.4, and m_n = 2 × 1.66 × 10⁻²⁴ g. We also assume that the ion density ρ_i is determined by the approximate relation (Elmegreen 1979; Nakano 1979)

$$\begin{equation} \rho _{\rm i}=m_{\rm i} K \left(\frac{\rho /m_{\rm n}}{10^5\, \mbox{cm$^{-3}$}}\right)^k, \end{equation} \tag{ 9 }$$

where we assume k = 0.5 throughout this paper, and take K = 3 × 10⁻³ cm⁻³ as a typical value. By using Equation (9), Equation (8) can be written as

$$\begin{equation} \tau _{\rm ni}=\gamma \rho ^{-1/2}, \end{equation} \tag{ 10 }$$

where γ ≃ 170.2 g^1/2 cm^−3/2 s using the above typical values. Then, Equation (3) becomes

$$\begin{equation} \frac{\partial \mbox{\boldmath $B$}}{\partial t} = \nabla \times (\mbox{\boldmath $v$} \times \mbox{\boldmath $B$}) + \nabla \times \left[\frac{\gamma }{c\rho ^{3/2}} (\mbox{\boldmath $j$} \times \mbox{\boldmath $B$}) \times \mbox{\boldmath $B$}\right]. \end{equation} \tag{ 11 }$$

As an initial condition, we assume hydrostatic equilibrium of a self-gravitating one-dimensional cloud along the z-direction (Kudoh & Basu 2003, 2006). The hydrostatic equilibrium is calculated from the equations

$$\begin{equation} \frac{dp}{dz}=\rho g_z, \end{equation} \tag{ 12 }$$

$$\begin{equation} \frac{dg_z}{dz}=-4\pi G \rho, \end{equation} \tag{ 13 }$$

$$\begin{equation} p=c_s^2 \rho, \end{equation} \tag{ 14 }$$

subject to the boundary conditions

$$\begin{equation} g_z(z=0)=0,\ \ \rho (z=0)=\rho _0,\ \ p(z=0)=\rho _0\, c_{s0}^2, \end{equation} \tag{ 15 }$$

where ρ₀ and c_s0 are the initial density and sound speed at z = 0. If the initial sound speed (temperature) is uniform throughout the region, we have the following analytic solution ρ_S found by Spitzer (1942):

$$\begin{equation} \rho _S(z)=\rho _0 \,\mbox{sech}^2(z/H_0), \end{equation} \tag{ 16 }$$

where

$$\begin{equation} H_0=\frac{c_{s0}}{\sqrt{2\pi G \rho _0}} \end{equation} \tag{ 17 }$$

is the scale height. However, an isothermal molecular cloud is usually surrounded by warm material, such as neutral hydrogen gas. Observations also show that the transition between molecular gas and surrounding gas is quite sharp (Blitz 1991). In order to simulate the situation, we assume the initial sound speed distribution to be

$$\begin{equation} c_{s}^2(z)=c_{s0}^2 + \frac{1}{2} (c_{sc}^2 - c_{s0}^2) \left[ 1+\tanh \left(\frac{|z|-z_c}{z_d}\right) \right], \end{equation} \tag{ 18 }$$

where we take c²_sc = 10 c²_s0, z_c = 2 H₀, and z_d = 0.1 H₀ throughout the paper. This equation gives a numerical model of the temperature transition of c²_sc/c²_s0 at z_c with the transition length of z_d. By using this sound speed distribution, we can solve Equations (12)–(14) numerically. The initial density distribution of the numerical solution shows that it is almost the same as Spitzer's solution for 0 ⩽ z ⩽ z_c.

We also assume that the initial magnetic field is uniform along the z-direction:

$$\begin{equation} B_z=B_0,\ \ B_x=B_y=0, \end{equation} \tag{ 19 }$$

where B₀ is constant.

In this equilibrium sheet-like gas, we input a random velocity fluctuation at each grid point:

$$\begin{equation} v_x=v_a R_m(x,y),\ \ v_y=v_a R_m(x,y),\ \ v_z=0.0, \end{equation} \tag{ 20 }$$

where R_m is a random number with spectrum v²_k ∝ kⁿ in Fourier space, where k = (k²_x + k²_y)^1/2, and the root mean square value is about 1. The R_m's for each of v_x and v_y are independent realizations. In this paper, we apply the two distinct spectra of v²_k ∝ k⁻⁴ and v²_k ∝ k⁰. The parameter v_a shows the strength of the velocity fluctuation. Models with the same spectrum use the same pair of realizations of R_m for generating the initial perturbations.

The constant G along with the two parameters c_s0 and ρ₀ that are associated with the initial state allow one to choose a set of three fundamental units for this problem. We choose these to be c_s0, H₀, and ρ₀ for velocity, length, and density, respectively. These yield a time unit t₀ ≡ H₀/c_s0. The initial magnetic field strength introduces one dimensionless free parameter,

$$\begin{equation} \beta _{0} \equiv \frac{8 \pi p_0}{B_0^2} = \frac{8 \pi \rho _0 c_{s0}^2}{B_0^2}, \end{equation} \tag{ 21 }$$

the initial ratio of gas to magnetic pressure at z = 0.

In the sheet-like equilibrium cloud with a vertical magnetic field, β₀ is related to the mass-to-flux ratio for Spitzer's self-gravitating cloud. The mass-to-flux ratio normalized to the critical value is

$$\begin{equation} \mu _S \equiv 2\pi G^{1/2} \frac{\Sigma _S}{B_0}, \end{equation} \tag{ 22 }$$

where

$$\begin{equation} \Sigma _S=\int _{-\infty }^{\infty } \rho _S dz = 2\rho _0 H_0 \end{equation} \tag{ 23 }$$

is the column density of Spitzer's self-gravitating cloud. Therefore,

$$\begin{equation} \beta _{0}=\mu _S^2. \end{equation} \tag{ 24 }$$

Although the initial cloud we used is not exactly the same as the Spitzer cloud, β₀ is a good indicator to whether or not the magnetic field can prevent gravitational instability in the flux-freezing case (Nakano & Nakamura 1978).

We also define a dimensionless parameter

$$\begin{equation} \alpha \equiv \gamma \sqrt{2 \pi G}, \end{equation} \tag{ 25 }$$

which shows the strength of the ambipolar diffusion. We note that α can also be equated to τ_ni,0/t₀, where τ_ni,0 is the initial value of the neutral–ion collision time calculated from Equation (10). By using the typical values in Section 2.1, we get α ≃ 0.11, which is taken as a fiducial parameter in this paper. We also vary α as a free parameter in Section 3.3.

Dimensional values of all quantities can be found through a choice of ρ₀ and c_s0. For example, for c_s0 = 0.2 km s⁻¹ and n₀ = ρ₀/m_n = 10⁴ cm⁻³, we get H₀ ≃ 0.05 pc, t₀ ≃ 2.5 × 10⁵ yr, and B₀ ≃ 20 μG if β₀ = 1.

In order to solve the equations numerically, we use the CIP method (Yabe & Aoki 1991; Yabe et al. 2001) for Equations (1), (2), and (7), and the method of characteristics-constrained transport (MOCCT; Stone & Norman 1992) for Equation (3), including an explicit integration of the ambipolar diffusion term. The combination of the CIP and MOCCT methods is summarized in Kudoh et al. (1999) and Ogata et al. (2004). It includes the CCUP method (Yabe & Wang 1991) for the calculation of gas pressure in order to get more numerically stable results. The numerical code used for this paper is based on that of Ogata et al. (2004).

In this paper, the ambipolar diffusion term is only included when the density is greater than a certain value, ρ_cr. We let ρ_cr = 0.3ρ₀ both for numerical convenience and due to the physical idea that the low density region is affected by external ultraviolet radiation so that the ionization fraction becomes large, i.e., τ_ni becomes small (Ciolek & Mouschovias 1995). Under this assumption, the upper atmosphere of the sheet-like cloud is not affected by ambipolar diffusion. This simple assumption helps to avoid very small time steps due to the low density region in order to maintain stability of the explicit numerical scheme.

We use a mirror-symmetric boundary condition at z = 0 and periodic boundaries in the x- and y-directions. At the upper boundary, at z = z_out = 4H₀, we also use a mirror-symmetric boundary except when we solve for the gravitational potential. This symmetric condition is just for numerical convenience. When the same boundary conditions were used by Kudoh et al. (2007) and Kudoh & Basu (2008) for three-dimensional simulations, the results were consistent with two-dimensional thin-disk simulations (e.g., Basu & Ciolek 2004; Basu et al. 2009a). Hence, we believe that the boundary conditions do not affect the result significantly. The Poisson Equation (5) is solved by the Green's function method to compute the gravitational kernels in the z-direction, along with a Fourier transform method in the x- and y-directions (Miyama et al. 1987b). This method of solving the Poisson equation allows us to find the gravitational potential of a vertically isolated cloud within |z| < z_out.

The computational region is |x|, |y| ⩽ 8πH₀ and 0 ⩽ z ⩽ 4 H₀. The number of grid points for each direction is (N_x, N_y, N_z) = (256, 256, 40). Since the most unstable wavelength for no magnetic field is about 4πH₀ (Miyama et al. 1987a), we have 64 grid points within this wavelength. We also have 10 grid points within the scale height of the initial cloud in the z-direction. The computational time of the fiducial model is about 70 hr of CPU time using four processors of the SX-9 in the National Astronomical Observatory of Japan. The maximum computational time, which occurs for the case of the highly subcritical cloud with small initial velocity fluctuations, is about 770 hr of CPU time. In the cases of the supercritical model, the computational times are about 30 minutes of CPU time.

We show the result of model V4 as a fiducial model, in which the initial turbulent velocity amplitude (v_a) is three times the sound speed (c_s0), its spectrum is v²_k ∝ k⁻⁴, the initial normalized mass-to-flux ratio is about 0.5 (i.e., β₀ = 0.25), and the dimensionless ambipolar diffusion coefficient has a typical value (α = 0.11). We suppose that these values are approximately typical in molecular clouds.

Figure 1 shows the time snapshot of the logarithmic density color map overlaid with the velocity vector field for model V4. The snapshot is at the end of our simulation, when the maximum density is greater than 100ρ₀. The upper panel shows the cross section at z = 0, and the lower panel shows it at y = 20.7 H₀. The maximum density is located at (x, y, z) = (−13.4 H₀, 20.7 H₀, 0). A collapsing core is located in the vicinity of the maximum density. Figure 2 shows the time snapshot of the logarithmic plasma β, the ratio of gas pressure to magnetic pressure, at the end of the simulation. It shows that β is greater than 1 around the core. It means that once the collapsing core is formed, the mass-to-flux ratio of the core is expected to be greater than 1 (β>1), i.e., supercritical. The other collapsing cores are also formed in the region where plasma β is greater than 1 (the left-down region of the upper panel in Figure 2). The left panel of Figure 3 shows the density, x-velocity, and plasma β along an x-axis cut at y = 20.7 H₀ and z = 0. The right panel shows the density, z-velocity, and plasma β along a z-axis cut at x = −13.4 H₀ and y = 20.7 H₀. In each panel, the velocity shows infall motion into the center of the core, though the x-velocity outside of the core (low density region) is affected by turbulent motions originating from the initial perturbation. The relative infall speed is subsonic and is about half of the sound speed at that time. The subsonic infall feature is consistent with previous studies of core formation in subcritical clouds (e.g., Ciolek & Basu 2000; Basu & Ciolek 2004).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 4 shows the time evolution of the maximum density at z = 0 and β at the location of the maximum density. The density goes up to ∼10 times greater than the initial density during the first compression of the cloud. Then, it rebounds and shows oscillations. Eventually, the dense region goes into runaway collapse around t = 18t₀, which corresponds to about 5 × 10⁶ yr. In the first compression, β also goes up to ∼0.6. It gradually goes up and becomes greater than 1 when the runaway collapse occurs. As we discussed in the previous section, β roughly equals the square of the normalized mass-to-flux ratio in the gravitational equilibrium state of gas with vertical magnetic field. The runaway collaping core is expected to be supercritical (β>1). Figures 2–4 demonstrate that β is a good indicator of the mass-to-flux ratio, though the gas is not exactly in gravitational equilibrium at the last stage of our simulations.

Zoom In Zoom Out Reset image size

Figure 5 shows a close-up view of the core in the vicinity of (x, y, z) = (−13.4 H₀, 20.7 H₀, 0). The iso-surface contour shows the logarithmic density and the lines represent the magnetic field. The core is located in a filamentary structure that was induced by the initial velocity fluctuation. The magnetic field lines show an hourglass-shaped structure because of the infall motion into the center of the core.

Zoom In Zoom Out Reset image size

Figure 6 shows a time snapshot of the logarithmic density overlaid with velocity vectors for model V0. The model V0 corresponds to the model for initial velocity amplitude (v_a) is 0.01 times sound speed (c_s0), but other parameters are the same as the fiducial model V4. The snapshot is at the end of our simulation, when the maximum density is greater than 100ρ₀. When the initial perturbed velocity is much smaller than the sound speed, the core formation time (t_core ∼ 136t₀, which corresponds to ∼3 × 10⁷ yr) is much longer than that (t_core ∼ 18.7t₀) of model V4. The maximum density is located at (x, y, z) = (7.6 H₀, − 11.9 H₀, 0). A collapsing core is formed in the vicinity of the maximum density. Figure 7 shows the time snapshot of the logarithmic plasma β at the end of the simulation for model V0. As we expected, the plasma β value is greater than 1 around the core, which means the core is supercritical. In Figures 6 and 7, core-like structures are gathered together on the side of x > 0. This reflects the initial perturbation of k⁻⁴ in which the larger scale contains greater energy in the velocity spectrum. The left panel of Figure 8 shows the density, x-velocity, and plasma β along an x-axis cut at y = 7.6 H₀ and z = 0. The right panel shows the density, z-velocity, and plasma β along a z-axis cut at x = −7.6 H₀ and y = −11.9 H₀. The velocity shows infall motion into the center of the core, though the center of the core shows a systematic speed in the positive x-direction of about 0.3 c_s0. The relative infall speed is also subsonic. Figure 9 shows the close-up view of the core in the vicinity of (x, y, z) = (−13.4 H₀, 20.7 H₀, 0). The iso-surface contour shows the logarithmic density and the lines represent the magnetic field. The core shows typical oblate-like structure. The magnetic field line shows a hourglass-shaped structure and the deformation of the field line is similar to that of model V4 in Figure 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 10 shows the time snapshot of the logarithmic density overlaid with velocity vectors at the end of the simulation for model K1. The model K1 is one for which the initial velocity spectrum is v²_k ∝ k⁰, but other parameters are the same as the fiducial model V4. The flat spectrum is possibly not consistent with observations, but it is used to compare with the result from the typical spectrum of v²_k ∝ k⁻⁴. When the initial perturbed velocity has a flat spectrum, the core formation time (t_core ∼ 67t₀) is longer than that (t_core ∼ 18.7t₀) of model V4. The maximum density is located at (x, y, z) = (−1.5 H₀, − 9.5 H₀, 0) and a collapsing core is formed in the vicinity of the maximum density. Figure 11 shows the time snapshot of the logarithmic plasma β at the end of the simulation for model K1. We note that β is greater than 1 around the core again. The left panel of Figure 12 shows the density, x-velocity, and plasma β along an x-axis cut at y = −9.5 H₀ and z = 0. The right panel shows the density, z-velocity, and plasma β along a z-axis cut at x = −1.5 H₀ and y = −9.5 H₀. The velocity shows infall motion into the center of the core, and the relative infall speed is subsonic again. Figure 13 shows a close-up view of the core in the vicinity of (x, y, z) = (−1.5 H₀, − 9.5 H₀, 0). The magnetic field lines also show an hourglass-shaped structure. Figures 3, 8, and 12 show the structure of collapsing cores are not so dependent on the initial velocity fluctuations, except that the outer region of the core in model V4 is strongly affected by the large-scale turbulent flow.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The main difference of each model is the core formation time. Figure 14 shows the time evolution of the maximum density for models V4, V0, and K1. In the case of model V4, which has a k⁻⁴ spectrum, the maximum density is strongly increased by the compression caused by the large-scale supersonic flow. In the case of model K1, the compression occurs in the initial stage, but it is weaker than that of model V4 because the small-scale perturbations contain more power than in model V4. In the case of model V0, with small-amplitude initial perturbations, the core formation time is about seven times longer than that of model V4.

Zoom In Zoom Out Reset image size

Figure 15 shows the core formation time as a function of the strength of the initial velocity perturbation. As the perturbed velocity amplitude increases, the core formation time decreases for each type of initial spectrum. In the cases of the initial k⁻⁴ spectrum (filled squares), the core formation times are systematically shorter than those of the initial k⁰ spectrum (open squares). In order to understand the physics of the core formation time, we plot the core formation time as a function of ρ_peak in Figure 16, where ρ_peak is defined as the value of the density peak during the first compression in the time evolution of the maximum density on z = 0. In the cases of the k⁻⁴ spectrum (filled squares), ρ_peak is greater than in the cases of the k⁰ spectrum (open squares), even when v_a is relatively small. Figure 16 shows that the core formation time is shorter when the density peak ρ_peak is greater, and indicates that it is nearly proportional to $1/\sqrt{\rho _{{\rm peak}}}$. The density dependence is similar to that derived in quasistatically contracting magnetically supported cores, assuming a force balance of the magnetic force and gravity in the core (e.g., Mouschovias & Ciolek 1999). In our simulation, the density is increased by the compression caused by large-scale nonlinear flows, and a dense region undergoes a rebound and oscillations, which means the core is not in an exact force balance. Nevertheless, it is interesting that the core formation time is nearly proportional to $1/\sqrt{\rho _{{\rm peak}}}$. We discuss this subject further in Section 4.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 17 shows the time evolution of the maximum density for models V4, A0, A1 and A2. The models A0, A1, and A2 correspond to the models for α = 0, α = 0.05, and α = 0.2, respectively, but their other parameters are the same as those of V4.

Zoom In Zoom Out Reset image size

As α decreases, the core formation time increases. When α = 0 (model A0), which means no ambipolar diffusion, a core was not formed during the calculation time, lasting until t ≃ 268t₀. Figure 18 shows the core formation time as a function of α. When the initial velocity perturbation is small, the core formation time shows nearly linear dependence on α⁻¹ (open triangles). When the initial perturbation is nonlinear (filled triangles), the dependence on α⁻¹ is a little steeper than a linear dependence. For model A2 (α = 0.2, v_a = 3c_s0), the core formation occurs very rapidly, without showing an oscillation of the maximum density. When α is large enough, the supercritical region can appear quickly after the first compression, resulting in a significantly reduced core formation time.

Zoom In Zoom Out Reset image size

In Figure 19, we show the time evolution of the total kinetic energies in the computational region for models V4, A0, A1, and A2. The oscillation of the energies comes from the strong compression and rebound of the clouds. The figure shows that the kinetic energies quickly decrease within a crossing time of the initial perturbed velocity across the computational region. The e-folding time of the energies depends slightly on α. As α increases, the e-folding time increases. In the case of α = 0 (model A0), about 10% of the initial kinetic energy remains at the end of the calculation (t ≃ 268t₀). Basu & Dapp (2010) found that a significant portion (∼1/2) of the initial kinetic energy remains in the flux-freezing case using a thin-disk approximation. Our result shows that the kinetic energy is more dissipative than in the thin-disk approximation. This is because of the different boundary condition in the vertical direction. We discuss this subject further in Section 4.

Zoom In Zoom Out Reset image size

Figure 20 shows the time evolution of the maximum densities for models V4, B0, B1, and B2, showing the difference between models with different β₀. The models B0, B1, and B2 correspond to the models for β = 0.04, β = 0.09, and β = 0.36, respectively, but their other parameters are the same as those of V4. As we have shown in Section 2, β₀ approximately equals the square of the initial normalized mass-to-flux ratio of the cloud. The figure shows that the core formation time is longer when β₀ is smaller. The figure also shows that the initial density enhancement is smaller when β₀ is smaller. Figure 21 shows the core formation time as a function of β₀. When β₀ is greater than 1, the initial cloud is supercritical. In the case of the supercritical clouds, the core formation occurs on a dynamical timescale. When the initial velocity is large, a core forms during the first compression. The core formation time of the supercritical clouds does not strongly depend on the initial mass-to-flux ratio (or β₀). When β₀ is smaller than 1, the initial cloud is subcritical. The core formation times of a subcritical cloud are generally longer than those of a supercritical cloud. In the case of a highly subcritical cloud, however, the core formation time does not strongly depend on the initial mass-to-flux ratio. The transition occurs near the critical cloud, β₀ ∼ 1. This tendency is consistent with that of the growth time from the linear analysis using the thin-disk approximation (Ciolek & Basu 2006) and with that of the ambipolar diffusion time in quasistatically contracting magnetically supported cores (Mouschovias & Ciolek 1999). It is interesting that the same tendency is achieved even when the initial velocity fluctuation is nonlinear (filled circles).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

One of the problems for core formation in subcritical clouds is that the timescale of the core formation is longer than that expected from observations (e.g., Jijina et al. 1999). In order to solve this problem, Zweibel (2002) and Fatuzzo & Adams (2002) pointed out that the ambipolar diffusion rate can be enhanced in a turbulent medium with a fluctuating magnetic field. In a different way, Li & Nakamura (2004) and Nakamura & Li (2005) found that the compression of the cloud by the large-scale turbulent flow shortened the timescale of core formation even in subcritical clouds. Kudoh & Basu (2008) and this paper followed the latter idea using three-dimensional simulations without a thin-disk approximation. Here, we discuss how the compression shortens the timescale of core formation in subcritical clouds.

Figure 16 indicates that the core formation time is nearly proportional to $1/\sqrt{\rho _{{\rm peak}}}$, where ρ_peak is the value of the density peak during the first compression in the time evolution of the maximum density at z = 0. The ambipolar diffusion time (τ_ad) is estimated from Equation (11) as $\tau _{\rm ad} \sim 4\pi (\sqrt{2\pi G}/\alpha) \rho ^{3/2}L^2/B^2$, where L is the gradient length scale introduced by the turbulent compression. Because the compression by the nonlinear flow is nearly one dimensional, the magnetic field scales roughly as B ∝ L⁻¹ within the flux-freezing approximation. The surface density Σ also scales as Σ ∝ L⁻¹. If the compression is rapid enough that vertical hydrostatic equilibrium cannot established, then ρ ∝ Σ ∝ L⁻¹ and τ_ad ∝ L^5/2 ∝ ρ^−5/2 (Elmegreen 2007; Kudoh & Basu 2008). This means that diffusion can occur quickly if the turbulent compression creates small values of L or large values of ρ. This would lead to a rapidly rising value of β in Figure 4 at t/t₀ ∼ 1. If diffusion is so effective during the first turbulent compression that a dense region becomes supercritical, then it will evolve directly into collapse. Alternatively, the stored magnetic energy of the compressed (and still subcritical) region may lead to re-expansion of the dense region. If re-expansion of the initial compression does occur, then there is enough time for the vertical structure to settle back to near-hydrostatic equilibrium. In this case, the density scales roughly like ρ ∝ Σ² ∝ L⁻² and τ_ad ∝ L ∝ ρ^−1/2. The ambipolar diffusion time is also estimated by Mouschovias & Ciolek (1999) as τ_ad ∝ ρ^−1/2 in quasistatically contracting magnetically supported cores, assuming a radial force balance between gravity and the magnetic Lorentz force. In our simulation, force balance is not exactly achieved, but the time average (of the oscillations) in the cores can be approximately in force balance. Since the first strong compression leads to rapid magnetic flux loss, a higher density region (than the initial value) is rapidly attained and then settles into an approximate force balance. Since the phase of continuing ambipolar diffusion in the oscillating high density region takes longer than the initial compression, the overall timescale of the core formation is approximately proportional to ρ^−1/2_peak. Here, instead of the time-average density, ρ_peak is used as a representative density of the compressed gas, since it is complicated to determine the time-average density of the moving, oscillating, and collapsing high density region. It is interesting that the relation attains even when the ρ_peak is used. We expect that the average density of the compressed overdense region would be nearly proportional to the peak density.

The relation means that the density dependence of τ_ad ∝ ρ^−1/2 is the same as that of the free fall time. Figure 16 shows that the actual time of the core formation is about 30 times longer than the free fall time of gas with ρ_peak when α = 0.11 and β₀ = 0.5. The value of the dimensionless coefficient (∼30) is expected to be inversely proportional to α (Section 3.3), and not to be strongly dependent on the mass-to-flux ratio except when it is nearly critical (Section 3.4). Even when the time-average density is used instead of the peak density, the value of the dimensionless coefficient would not be so different because of the weak dependence of ρ (∝ρ^−1/2).

In some cases, collapse may occur during the first compression itself or soon thereafter, and the approximate scaling of the core formation time with ρ^−1/2_peak will not hold. Whether or not this occurs depends not only on the strength of the turbulent compression but also on the initial mass-to-flux ratio of the cloud (related to β₀) and the ambipolar diffusion coefficient (α). In model B2 (a subcritical model that is closest to critical, with normalized initial mass-to-flux ratio ≃0.6), a re-expansion does occur but the collapse starts relatively quickly in the third oscillation (Figure 20). In model A2 (α = 0.2), which has the poorest neutral–ion coupling of all models, the collapse starts almost in the first compression (Figure 17).

An overall conclusion is that the core formation occurs more rapidly than it would in the initial state due to the elevated value of ρ in the compressed but oscillating region. Since the core formation time is approximately proportional to ρ^−1/2, 10–100 times density enhancement is needed to get 3–10 times shorter core formation time than that of the standard core formation model in subcritical clouds (e.g., Jijina et al. 1999). The observed non-thermal velocity (3–10 times sound speed) is eligible to make such enhancement by compression in the isothermal clouds, if its scale is larger than the Jeans length.

Even though the core formation time is shortened by the nonlinear flows, the density, velocity, and magnetic field structure of a core do not strongly depend on the initial strength of the velocity fluctuation (e.g., Figures 3, 8, and 12). The infall velocities of the cores are subsonic and the magnetic field lines show weak hourglass shapes. This result may be consistent with the result that the core formation time is proportional to ρ^−1/2_peak. Even when core formation is initiated by the initial turbulence, the core properties can be similar to the quasistatically contracting magnetically supported cores discussed by Mouschovias & Ciolek (1999). The initial turbulence accelerates core formation, but it eventually dissipates in the dense region when the collapsing core is formed. Subsonic infall motions were found by Lee et al. (2001) in an observational survey of starless cores. They found that the typical infall radii are 0.06–0.14 pc and that the infall speed lies in the range of 0.05–0.09 km s⁻¹. These values are consistent with our results. The subsonic infall is the common feature of core formation in subcritical clouds, which has been pointed out by one-dimensional axisymmetric (Ciolek & Basu 2000) and two-dimensional thin-disk (Basu & Ciolek 2004; Basu et al. 2009a, 2009b) models. Because of the subcritical infall, the hourglass shapes are a little weaker than that of an initially supercritical cloud that shows transsonic infall (Kudoh et al. 2007). More quantitative analysis of the field morphological difference and the infall speeds may distinguish the mass-to-flux ratio in a molecular cloud.

In Figure 19, we showed the time evolution of the total turbulent kinetic energy in isothermal subcritical clouds. It shows that the kinetic energy dissipates efficiently even in the flux-freezing limit (α = 0). This is different from the result of Basu & Dapp (2010). They found that the turbulent kinetic energy persists in an isothermal subcritical cloud in the flux-freezing limit. They showed that the energy dissipation occurs effectively only in the stage of the first nonlinear compression. After that, the turbulent energy persists, and oscillates about an average value. The average value is estimated to be about half of the initial turbulent kinetic energy for the same parameters as those of our fiducial model. On the other hand, our simulation shows that energy dissipation occurs effectively even in the flux-freezing limit; energy dissipation occurs almost exponentially during several oscillations, although about 12% of the initial kinetic energy remains at the final stage. The difference is likely due to the different vertical boundary conditions on the molecular cloud. Basu & Dapp (2010) use a thin-disk approximation of the cloud, and adopt the current-free condition of the magnetic field outside of the disk, assuming the outside density is negligibly small. We study the problem by three-dimensional MHD and without a thin-disk approximation, but assume that the molecular cloud is surrounded by warm gas whose density is 10 times smaller than the molecular cloud. The computational region in our model is 0 < z < 4 H₀, and so does not cover a very large dynamic range of densities. Therefore, future three-dimensional simulations with a larger computational region along the z-direction and including low-density gas outside of the disk may make for more revealing comparisons with the result of Basu & Dapp (2010).

In addition to the parameter study in Table 1, we performed simulations of models with the same parameters as models V4 and V1, but different random realizations of the initial velocity fluctuations. The core formation time for different realizations varies by about 25% in these samples. The overall evolution and the result of the accelerated core formation by nonlinear flows were not altered by the different initial random realizations. We performed simulations with different spatial resolutions for the model V4. In the case of the highest resolution of (N_x, N_y, N_z) = (512, 512, 40), the core formation time is t_core = 14.2t₀. There is a slight tendency that the core formation time becomes a little shorter for the high-resolution cases (Kudoh & Basu 2009). However, we should note that a realization of the random perturbations for the initial velocity fluctuation is also not the same for the different resolutions, because the random perturbations are input on all scales down to the grid scale. At the least, we can say that the result of the accelerated core formation by nonlinear flows is not altered significantly by spatial resolution. We also performed a simulation with a different boundary condition of magnetic fields for the model V4. Instead of the symmetrical boundary at z = z_out, we used the vertical field condition for the magnetic field, i.e., B_x = B_y = 0 for z ⩾ z_out. The overall evolution and the result of the accelerated core formation were not also altered by this. The core formation time for the different boundary varies by about 0.2% in this sample.