RECURRENT PLANET FORMATION AND INTERMITTENT PROTOSTELLAR OUTFLOWS INDUCED BY EPISODIC MASS ACCRETION

To study the formation and evolution of a circumstellar disk in a magnetized molecular cloud core, we solve the three-dimensional resistive MHD equations, including self-gravity:

$$\begin{eqnarray} & {\displaystyle \frac{\partial \rho }{\partial t}} + \nabla \cdot (\rho \mbox{\boldmath $v$}) = 0, & \end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray} &&\rho {\displaystyle \frac{\partial \mbox{\boldmath $v$}}{\partial t}} + \rho (\mbox{\boldmath $v$} \cdot \nabla)\mbox{\boldmath $v$} = - \nabla P - {\displaystyle \frac{1}{4 \pi }} \mbox{\boldmath $B$} \times (\nabla \times \mbox{\boldmath $B$}) - \rho \nabla \phi,\quad \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} & {\displaystyle \frac{\partial \mbox{\boldmath $B$}}{\partial t}} = \nabla \times (\mbox{\boldmath $v$} \times \mbox{\boldmath $B$}) + \eta _{\rm OD} \nabla ^2 \mbox{\boldmath $B$}, & \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} & \nabla ^2 \phi = 4 \pi G \rho, & \end{eqnarray} \tag{ 4 }$$

where ρ, $\mbox{\boldmath $v$}$, P, $\mbox{\boldmath $B$}$, η_OD, and ϕ denote density, velocity, pressure, magnetic flux density, Ohmic resistivity, and gravitational potential, respectively. In addition, we adopted the hyperbolic divergence cleaning method of Dedner et al. (2002) to obtain a divergence-free magnetic field ($\mbox{\boldmath $\nabla $} \cdot \mbox{\boldmath $B$}=0$). With this method, no magnetic monopoles appear throughout the calculation (see also Machida et al. 2005b; Matsumoto 2007).

To mimic the temperature evolution calculated by Masunaga & Inutsuka (2000), we adopt the piecewise polytropic equation of state (see Vorobyov & Basu 2006; Machida et al. 2007) as

$$\begin{equation} P = c_{s,0}^2\, \rho \left[ 1+ \left({\displaystyle \frac{\rho }{\rho _c}} \right)^{2/5} \right], \end{equation} \tag{ 5 }$$

where c_s,0 = 190 m s⁻¹ and ρ_c = 3.84 × 10⁻¹⁴ g  cm⁻³ (n_c = 10¹⁰  cm⁻³). With Equation (5), the gas behaves isothermally for n ≲ 10¹⁰  cm⁻³ and adiabatically for n ≳ 10¹⁰  cm⁻³. For a realistic evolution of the magnetic field in the circumstellar disk, we adopt the resistivity (η_OD) as the fiducial value in Machida et al. (2007), in which the Ohmic dissipation becomes effective for 10¹¹  cm⁻³ ≲ n ≲ 10¹⁵  cm⁻³ (for details, see Equations (9) and (10), and Figure 1 of Machida et al. 2007).

As the initial state, we take a spherical cloud core with a critical Bonnor–Ebert (BE) density profile, in which a uniform density (ρ_amb ≃ 0.07ρ_c) is adopted outside the sphere (r>R_c). For the BE density profile, we adopt the central density as n_c = 10⁶  cm⁻³ and isothermal temperature as T = 10 K. For these parameters, the critical BE radius is R_c = 4.7 × 10³ AU. The gravitational force is ignored outside the host cloud (r>R_c) to mimic a stationary interstellar medium. In addition, we prohibit the gas inflow at r = R_c from suppressing mass input from the interstellar medium into the gravitationally collapsing cloud core. Thus, only the gas inside r < R_c collapses to form the circumstellar disk and protostar. Note that although the protostellar outflow driven by the circumstellar disk propagates into the region of r>R_c and disturbs the interstellar medium over time, we can safely calculate the mass accretion process onto the circumstellar disk because the outflow propagating into the interstellar medium does not affect the inner cloud core region (r < R_c). In this study, we call this initial spherical cloud with BE density profile "the cloud core" that has a radius of R_c and mass of M_cl.

Since the critical BE sphere is in an equilibrium state, we increase the density by a factor of f to promote the contraction, where f is the density enhancement factor that represents the stability of the initial cloud core. An initial cloud core with larger f is more unstable against gravity. In general, the stability of the cloud core is represented by a parameter α₀ (≡E_th/E_grav), which is the ratio of thermal (E_th) to gravitational (E_grav) energy. As shown in Matsumoto & Hanawa (2003), when the BE density profile is adopted, the density enhancement factor is related to the parameter α₀ as

$$\begin{equation} \alpha _0 = {\displaystyle \frac{0.84}{f}}. \end{equation} \tag{ 6 }$$

We constructed two models with different f (=1.68 and 2.8). In this paper, we call the model having f = 1.68 (α = 0.5) "model A05" and the model having f = 2.8 (α = 0.3) "model A03." At the initial state, model A03 is more unstable than model A05.

In both models, the initial cloud core rotates rigidly around the z-axis in the region of r < R_c, while the uniform magnetic field parallel to the z-axis (or rotation axis) is adopted in the whole computational domain. In addition, the angular velocity decreasing in proportion to exp(−r²) outside the host cloud (r>R_c) is also adopted. Both models have the same angular velocity with Ω₀ = 1.1 × 10⁻¹³ s⁻¹ and the same magnetic field, B₀ = 37 μG. The mass inside r < R_c for each model is 0.8 M_☉ (model A05) and 1.3 M_☉ (model A03). Model names and initial values are summarized in Table 1. The initially magnetized, rotating cloud core is also specified by the non-dimensional parameters β₀ and γ₀, where β₀ (≡E_rot/|E_grav|) and γ₀ (≡E_mag/|E_grav|) is the ratio of the rotational and magnetic energy to the gravitational energy inside the initial cloud core. Model A03 has β₀ = 2 × 10⁻³ and γ₀ = 0.02, while model A05 has β₀ = 3 × 10⁻³ and γ₀ = 0.05. In addition, the mass-to-flux ratio M/Φ is also used to specify the initial cloud, where M and Φ are the mass and magnetic flux of the initial cloud core. There exists a critical value of M/Φ below which a cloud is supported against the gravity by the magnetic field. For a cloud with uniform density, Mouschovias & Spitzer (1976) derived a critical mass-to-flux ratio

$$\begin{equation} \left({\displaystyle \frac{M}{\Phi }} \right)_{\rm cri} = {\displaystyle \frac{\zeta }{3\pi }} \left({\displaystyle \frac{5}{G}} \right)^{1/2}, \end{equation} \tag{ 7 }$$

where the constant ζ = 0.53. We define the mass-to-flux ratio normalized by the critical value as λ. In our setting, model A03 has λ = 8.9, while model A05 has λ = 5.5.

Table 1. Model Parameters and Calculation Results

Model	α0	f	Ω0 (s−1)	Bini (μG)	Mcl ($\thinspace M_\odot$)	Rc (AU)	λ	Mps ($\thinspace M_\odot$)	Mdisk ($\thinspace M_\odot$)	Npl
A03	0.3	2.8	1.1 × 10−13	37	1.3	4700	8.9	0.59	0.15	16
A05	0.5	1.7	1.1 × 10−13	37	0.8	4700	5.5	0.15	0.03	0

Download table as:  ASCIITypeset image

We add a m = 2-mode non-axisymmetric density perturbation to the initial cloud core. Then, the density profile of the cloud core is described as

$$\begin{eqnarray} \rho (r) = \left\lbrace \begin{array}{@{}ll}\rho _{\rm BE}(r) \, (1+\delta _\rho)\,f & \mbox{for} \; \; r < R_{c}, \\ \rho _{\rm BE}(R_c)\, (1+\delta _\rho)\,f & \mbox{for}\; \; r \ge R_{c}, \\ \end{array} \right. \end{eqnarray} \tag{ 8 }$$

where ρ_BE(r) is the density distribution of the critical BE sphere, and δ_ρ is the axisymmetric density perturbation. For the m = 2-mode, we chose

$$\begin{equation} \delta _\rho = A_{\phi } (r/R_{\rm c})^2\, {\rm cos}\, 2\phi, \end{equation} \tag{ 9 }$$

where A_ϕ (= 0.01) represents the amplitude of the perturbation. The radial dependence is chosen so that the density perturbation remains regular at the origin (r = 0) at one timestep after the initial stage. This perturbation ensures that the center of gravity is always located at the origin.

In the collapsing cloud core, we assume the protostar formation to occur when the number density exceeds n>10¹³  cm⁻³ at the center. To model the protostar, we adopt a sink around the center of the computational domain. In the region r < r_sink = 1 AU, gas having a number density of n>10¹³  cm⁻³ is removed from the computational domain and added to the protostar as gravity in each timestep (for details, see Machida et al. 2009a, 2010a). This treatment of the sink makes it possible to calculate the evolution of the collapsing cloud core and circumstellar disk for a longer duration. In addition, inside the sink, the magnetic flux is removed by Ohmic dissipation, because such a region has a magnetic Reynolds Re number exceeding unity Re>1 (for details, see Machida et al. 2007).

To calculate on a large spatial scale, the nested grid method is adopted (for details, see Machida et al. 2005a, 2005b). Each level of a rectangular grid has the same number of cells (128 × 128 × 32). The calculation is first performed with five grid levels (l = 1–5). The box size of the coarsest grid l = 1 is chosen to be 2⁵R_c. Thus, the grid of l = 1 has a box size of ∼1.5 × 10⁵ AU. A new finer grid is generated before the Jeans condition is violated (Truelove et al. 1997). The maximum level of grids is l_max = 12 which has a box size of 74 AU and a cell width of 0.58 AU.

We adopted a fixed boundary condition on the outermost grid boundary. The uniform density (ρ_amb) and magnetic field (B_x = 0, B_y = 0, B_z = B₀) and zero fluid velocity (v_x = v_y = v_z = 0) are imposed on the l = 1 grid boundary at each timestep. Such a boundary condition hardly affects the evolution of the cloud core, because the gravitationally collapsing cloud core with a radius of 4.7 × 10³ AU (r = R_c) is embedded in a large simulation box with a size of ∼1.5 × 10⁵ AU, inside which the static interstellar medium is fulfilled in the region of r>R_c. In our setting, the Alfv$\acute{\rm e}$n speed in the interstellar medium (r>r_BE) for model A03 is v_A = 0.15 km s⁻¹. Thus, for the Alfv$\acute{\rm e}$n wave to reach the computational boundary from the center of the cloud core, it will take 2.4 × 10⁶ yr. Since we stopped our calculation ∼10⁴ yr after the calculation began, the Alfv$\acute{\rm e}$n wave generated at the center of the cloud core (or the computational boundary) never reaches the computational boundary (or the center of the cloud core).

As described in Section 2, we constructed two models (models A03 and A05) with different initial stabilities (i.e., different α₀). In the two models, the protostar forms t ≃ 2.50 × 10⁴ yr (model A03) and t ≃ 5.48 × 10⁴ yr (model A05) after the calculation begins. Thus, in the collapsing cloud core, the protostar for model A03 forms in a shorter duration than that for model A05. This is because the initial cloud core for model A03 is more unstable than that for model A05 and begins to collapse according to the self-similar solution (Larson 1969) in a shorter time after the calculation begins (Machida et al. 2008b). In both models, we calculated the evolution of the circumstellar disk for ∼10⁴ yr after protostar formation.

First, we describe the evolution of the circumstellar disk for model A03 that shows the recurrent fragmentation and formation of substellar-mass object in the main accretion phase. Figure 1 shows the evolution of the circumstellar disk after protostar formation for model A03, in which the density distributions around the protostar on the equatorial plane are plotted. The box sizes in the upper panels are ∼80 AU, while those in the middle and lower panels are ∼320 AU. Figure 1(a) shows the structure of the circumstellar disk t_c ≃ 450 yr after protostar formation, in which the red region (n ≳ 10¹⁰ − 10¹¹  cm⁻³) corresponds to the circumstellar disk. Note that we describe the elapsed time after protostar formation with the notation "t_c," which differs from the elapsed time "t" after the calculation begins (or the cloud begins to collapse). The size of the circumstellar disk increases with time and extends up to r∼ 200 AU by the end of the calculation. In the circumstellar disk, fragmentation occurs due to gravitational instability t_c ≃ 795 yr after protostar formation (for details, see Inutsuka et al. 2010 and Machida et al. 2010a). Two ambiguous clumps appear in Figures 1(b) and (c), while two clear fragments are seen in Figure 1(d). Fragmentation occurs ∼6.8 AU away from the protostar. Then, the fragments orbit around the protostar with an orbital separation of ∼2–38 AU (Figures 1(c)–(e)) and finally fall onto the protostar at t_c ≃ 3657 yr (Figures 1(e) and (f)). In this section, we call the fragment that appeared in the circumstellar disk the planet. Note that in reality, since some fragments that appeared in this calculation exceeded the deuterium-burning limit (about 13 Jupiter mass), they are not planets in the complete sense. Note also that we define, for convenience, the objects born in the circumstellar disk as planets to characterize them as a whole: we discuss their mass in Section 4.2.

Zoom In Zoom Out Reset image size

Figure 2 shows the structure of the circumstellar disk and protostellar outflow at almost the same epoch as Figure 1(c) in three dimensions. The outflow driven by the first core in the collapsing cloud core has been studied by many authors (e.g., Matsumoto & Tomisaka 2004; Machida et al. 2004, 2005b; Hennebelle & Fromang 2008). Even in the main accretion phase (i.e., even after protostar formation), the circumstellar disk that originates from the first core continues to drive the outflow (Tomisaka 2002; Banerjee & Pudritz 2006; Machida et al. 2007, 2009a). In this model, the outflow driven by the circumstellar disk extends up to ∼516 AU by the end of the calculation. As seen in Figure 2, the outflow is driven only by the outer disk region, while the outflow is not driven by the inner region of the circumstellar disk where fragmentation occurs and planets appear. This is because the inner disk region has a density of n ≳ 10¹¹–10¹²  cm⁻³ and the magnetic field is dissipated by Ohmic dissipation. As shown in Nakano et al. (2002) and Machida et al. (2007), when the number density exceeds n ≳ 10¹¹–10¹²  cm⁻³, Ohmic dissipation becomes effective owing to the extremely low ionization degree. Note that the ambipolar diffusion can be effective in the range of n < 10¹¹–10¹²  cm⁻³ (see Section 4.3). Fragmentation is suppressed by the magnetic field because magnetic effects such as magnetic braking and outflows effectively transfer the angular momentum that promotes fragmentation (Machida et al. 2005a, 2008b; Hennebelle & Teyssier 2008). Thus, a weaker field promotes fragmentation in the circumstellar disk. As a result, the inner disk region has a considerable weak magnetic field and fragmentation occurs without outflow, while the outer disk region has a strong field and no fragmentation occurs with outflow.

Zoom In Zoom Out Reset image size

In Figure 2, yellow lines are magnetic field lines. In the figure, only magnetic field lines that threaded the planets are plotted to stress the planet's spin or orbital motion. The magnetic field lines just above each planet are strongly twisted by the spin motion of the planet. In midair (or on a large scale), the magnetic field lines originating from each planet are twisted together, reflecting the orbital motion of the planets. Although we only plotted magnetic field lines near the planet to understand the planet's motion, the magnetic field (or Lorentz force) in this region (or inner disk region) is extremely weak (β_p ≫ 10, where β_p is the plasma beta and different from initial ratio of the rotational to the gravitational energy, β₀) because Ohmic dissipation is effective. Thus, near the planets, the magnetic field lines are passively moved, and they hardly affect the dynamical evolution of both the planet and circumstellar disk. However, the planets perturb the outer disk region that is connected to the magnetic field lines driving the outflow (i.e., the magnetic field lines in the outer disk region where Ohmic dissipation is not effective). As a result, a highly time-variable outflow appears in this model. We will discuss an intermittent outflow in Section 3.7.

Figure 1(e) shows the structure of the circumstellar disk just before the planets fall onto the protostar, while Figure 1(f) shows the disk just after the fall. After the first-generation planets disappear, fragmentation occurs again in the circumstellar disk and the second-generation planets appear t_c ≃ 3990 yr after protostar formation as shown in Figure 1(g). However, these planets also fall onto the protostar ∼680 yr after their formation. The formation and falling of the planets are repeated several times in the circumstellar disk. By the end of the calculation, 18 planets appeared and 16 planets fell onto the protostar. Figure 3 shows eighth- and ninth-generation planets. The eighth-generation planets fell onto the protostar, while the ninth-generation planets survived until the end of the calculation. The orbital radius r_i, time t_i at planet formation, the planet's falling epoch t_e, lifetime t_e − t_i, and mass at its formation M_pl are summarized in Table 2.

Zoom In Zoom Out Reset image size

Figure 4 shows the orbital evolution of planets against time after protostar formation t_c. When multiple planets appear as seen in Figures 1 and 3, only the orbital radius of the most massive planet is plotted. The figure indicates that the first-generation planets form ∼6.8 AU from the protostar at t_c ≃ 795 yr, and shrink their orbit and fall onto the protostar after ∼10 orbital rotations. The second-generation planets appear at ∼46 AU and fall onto the protostar after ∼2–3 orbital rotations. This formation and falling of planets is repeated until the ninth-generation planet is formed.

Zoom In Zoom Out Reset image size

Figure 5 shows the trajectory of the planets for each generation. The figure indicates that planets shrink their orbits with a certain amount of eccentricity. Figures 4 and 5 show that planets are born ∼6–50 AU from the protostar and orbit for ∼500–1000 yr (or 2–10 orbital rotations) around the protostar. In addition, these figures show that the first-generation planets are formed at a relatively small orbital radius (∼6.5 AU) and orbit for ∼10 times around the protostar, while other generations of planets are formed at a larger orbital radii (∼20–50 AU) and orbit only for 2–3 times. This is because the circumstellar disk increases its size and mass with time. As shown in Figure 1, the circumstellar disk has a size of ∼20–30 AU when the first-generation planet appears, while the disk size reaches 160 AU just before the second-generation planet appears. As a result, the magnetically inactive zone also expands as the disk mass and radius increase, and fragmentation and planet formation can occur in such an outer disk region. The mass evolution of the circumstellar disk is shown in Section 3.4, and the fragmentation condition is discussed in Section 3.6.

Zoom In Zoom Out Reset image size

Unlike model A03, model A05 shows no fragmentation in the circumstellar disk by the end of the calculation. Figure 6 shows the density distribution on the y = 0 (left) and z = 0 (right) planes t_c = 3781 yr after protostar formation for model A05. As in model A03, model A05 also shows the outflow driven by the circumstellar disk. The outflow appears ∼700 yr before protostar formation and reaches ∼4400 AU at the end of the calculation. In the right panel, the ring-like density gap is seen in the range of 40 AU ≲ r ≲ 80 AU. In this range, the gas is largely blown away from the circumstellar disk by the outflow.

Zoom In Zoom Out Reset image size

In the right panel, a non-axisymmetric structure is seen in the proximity of the protostar (r ≲ 40 AU), while nearly axisymmetric structures are seen in the outer disk region (r ≳ 40 AU). The magnetic field is well coupled with neutrals in the outer disk region because the gas density is lower than n ≲ 10¹¹  cm⁻³ where Ohmic dissipation is not effective. Thus, the angular momentum is effectively transferred by magnetic braking in the outer disk region, and the gas in such a region can fall into the inner disk region. On the other hand, the magnetic field is significantly dissipated by the Ohmic dissipation in the inner disk region (r ≲ 40 AU), because the density in such a region exceeds n ≳ 10¹¹  cm⁻³. Thus, in the inner disk region, angular momentum transfer by the magnetic field is not effective, and the excess angular momentum and accumulated mass cause gravitational instability to induce the non-axisymmetric structure (Walch et al. 2009a; Machida et al. 2010a). However, since the increase rate of the disk mass for model A03 is considerably smaller than that for model A05 (see Section 3.4), the circumstellar disk for model A05 is not sufficiently massive to induce fragmentation (see Section 3.6). For model A05, fragmentation and formation of a substellar-mass object do not occur until the end of the calculation. As a result, since the circumstellar disk is not disturbed by the fragments, the protostellar and circumstellar-disk mass moderately increase, unlike in the case of model A03 (see Section 3.4).

Figure 7 shows the evolution of protostellar mass against time after protostar formation. The protostellar mass M_ps is derived as the mass falling into the sink. At its formation, the protostar has a mass of $M_{\rm ps} \sim 10^{-3}\thinspace M_\odot$, as expected in spherically symmetric calculations (Larson 1969; Masunaga & Inutsuka 2000). In the main accretion phase, protostars increase their mass by gas accretion and reach M_ps = 0.59 M_☉ (model A03) and 0.15 M_☉ (model A05), respectively, at the end of the calculation. In Figure 7, the step-like increase of the protostellar mass for model A03 is caused by the falling of the planet. Since the planet and its envelope have a total mass of ∼6–27 M_Jup (see Table 2), the protostar drastically increases its mass at the moment of falling. In both models, the protostars have almost the same mass for t_c ≲ 3000 yr, while the protostar for model A03 is more massive than that for model A05 for t_c ≳ 3000 yr. The difference in protostellar masses is caused by the difference in the rate of mass accretion onto the protostar. The different mass accretion rates are due to different efficiencies of the angular momentum transfer in the circumstellar disk. The higher rate of mass accretion onto the protostar is realized in a more massive circumstellar disk, because a massive disk tends to show gravitational instability that effectively transfers the angular momentum outward, and increases the rate of mass accretion onto the protostar.

Zoom In Zoom Out Reset image size

Figure 8 shows the time evolution of the circumstellar-disk mass that is estimated in the same manner as in Machida et al. (2010a). The figure indicates that the disk mass for model A03 is significantly greater than that for model A05 at the epoch just after protostar formation. The circumstellar disk for model A05 gradually increases its mass until the end of the calculation and reaches M_disk ∼ 0.042 M_☉ for ∼10⁴ yr. The average increase rate of disk mass for model A05 is $\dot{M}_{\rm disk}\sim 4.2\times 10^{-6}\,M_\odot \,{\rm yr}^{-1}$. On the other hand, the circumstellar disk mass for model A03 reaches M_disk ∼ 0.15 M_☉ for ∼10⁴ yr with an average increase rate of $\dot{M}_{\rm disk}\sim 1.5\times 10^{-5}\,M_\odot \,{\rm yr}^{-1}$. Thus, considering that the initial stability of the cloud core is closely related to the mass accretion rate (see Section 4.2), the increase rate of the disk mass for model A03 is about one order of magnitude higher than that for model A05. Note that it is difficult to accurately estimate the total mass of accreting gas onto the circumstellar disk because a fraction of the circumstellar-disk mass falls onto the protostar with time. For model A03, the disk mass shows a high time variability that is caused by the time variable mass accretion onto the protostar (see Section 3.5); the spike-like decrease of the disk mass for model A03 corresponds to the epoch at which the planet falls onto the protostar.

Zoom In Zoom Out Reset image size

Since the gas accretes onto the protostar mainly through the circumstellar disk in the main accretion phase, the mass accretion rate depends on disk properties such as disk mass, spiral structure, and existence of fragments. Figure 9 shows the mass accretion rates for ∼10⁴ yr after protostar formation. We estimated the mass falling into the sink every year and defined it as the rate of mass accretion onto the protostar, $\dot{M}_{\rm ps}$. Thus, the time resolution of the accretion rate in Figure 9 is one year, while the timestep of calculations is Δt ≃ 0.001–0.01 yr. Therefore, we integrated the falling mass every ∼100–1000 timesteps.

Zoom In Zoom Out Reset image size

In Figure 9, the mass accretion rate for model A05 gradually decreases, oscillating between $\dot{M}\simeq 10^{-4}$ and $\dot{M} \sim 10^{-6}\,M_\odot \,{\rm yr}^{-1}$, while model A03 shows intermittent gas accretion with multiple peaks. For model A05, it is considered that the spiral structure in the inner disk region induces the non-steady accretion flow. The spiral structure can transfer the angular momentum outward and move the gas inward. The range of the accretion rate in this model almost corresponds to the theoretical prediction that is expressed as

$$\begin{equation} \dot{M}=f_{\rm m}\, {\displaystyle \frac{c_s^3}{G}} \simeq 1.9\times 10^{-6}\, f_{\rm m} \, M_\odot \,{\rm yr}^{-1}, \end{equation} \tag{ 10 }$$

where f_m is a numerical factor (e.g., f_m = 0.975 for Shu 1977, f_m = 46.9 for Hunter 1977), and c_s = 200 m s⁻¹ is adopted. The gas falls onto the circumstellar disk first, and the rate of gas accretion onto the circumstellar disk is ∼4 × 10⁻⁶ M_☉ yr⁻¹ on average, as shown in Figure 8. In summary, the rate of mass accretion onto the protostar almost corresponds to the mass accretion onto the circumstellar disk for model A05.

On the other hand, model A03 shows intermittent mass accretion caused by the large-scale perturbation of the disk and the falling of planets. As shown in Figure 4, the first-generation planets form at r ∼ 6.8 AU and orbit in the range of 2–40 AU until they fall onto the planet. They disturb the circumstellar disk and transfer the angular momentum outward. Thus, the mass accretion rate is linked to the orbital motion of the planet especially for t_c ≲ 3500 yr, i.e., before the first-generation planets fall onto the protostar. After the first-generation planets disappear, planets of subsequent generations (second to ninth) appear at r ≃ 20–50 AU from the protostar. They orbit in the range of 5–50 AU and approach the protostar within r < 5 AU only before falling onto the protostar. Thus, just before falling, these planets can disturb the inner disk region and promote mass accretion onto the protostar. The peak mass accretion rate for t_c ≳ 3500 yr corresponds to the epoch in which a planet approaches and falls onto the protostar. In addition, since not only the gas in the circumstellar disk but also the planet itself falls onto the protostar, the mass accretion rate at its peak for a model with planets is much larger than the theoretical prediction. In Figure 9, there are eight peaks of the mass accretion rate for t_c ≳ 3500 yr, corresponding to epochs of falling of planets. In such an epoch, the protostar has a very high mass accretion rate of $\dot{M}_{\rm ps} \gtrsim 10^{-3}\,\,\hbox{ to }10^{-2}\,M_\odot \,{\rm yr}^{-1}$, because a planet with mass of ${\sim} 10^{-3}\hbox{ to }10^{-2}\, M_\odot$ falls onto the protostar in a short duration. Note that this high accretion rate may be smoothed over time if calculated with a higher spatial resolution, while it is considered that the total mass falling onto the protostar does not change.

Figure 10 shows the radial distribution of the surface density σ (Figure 10(a)), Toomre Q-parameter (Figure 10(b); Toomre 1964), plasma β_p (Figure 10(c)), and ratio of the azimuthal to Kepler velocity v_ϕ/v_kep (Figure 10(c)). These quantities are azimuthally averaged. The azimuthal-averaged surface density σ [≡σ(r)] is derived as

$$\begin{equation} \sigma (r) = {\displaystyle \frac{dr \int _0^{2\pi } \sigma _{\rm s} (r, \theta)\, r_{\rm c}\, d\theta }{dS}}, \end{equation} \tag{ 11 }$$

where

$$\begin{equation} \sigma _{\rm s} (r, \theta) = \int ^{z < z_{\rm cri}} \rho (r, \theta, z)\, dz \end{equation} \tag{ 12 }$$

and

$$\begin{equation} {\it dS} = dr \int _0^{2\pi } r_{\rm c}\, d\theta, \end{equation} \tag{ 13 }$$

where r_c is the cylindrical radius and dr = ds(l) (the cell width of lth-level grid) is adopted. To determine z_cri in Equation (12), we estimated the disk height by eye and z_cri = 10 AU is adopted. We confirmed that the surface density σ hardly depends on z_cri in the range of z_cri>5 AU. In Figures 10(b) and (c), the Q parameter and plasma β_p are integrated along the z-axis and azimuthally averaged. To estimate the Q parameter, first, we integrated the mass-weighted sound speed c_s and epicyclic frequency κ along the z-axis as

$$\begin{equation} c_{\rm s,s} (r, \theta) = {\displaystyle \frac{\int ^{z < z_{\rm cri}} c_s(r, \theta, z)\, \rho (r, \theta, z)\,dz}{\sigma _{\rm s}(r, \theta)}} \end{equation} \tag{ 14 }$$

and

$$\begin{equation} \kappa _{\rm s} (r, \theta) = {\displaystyle \frac{\int ^{z < z_{\rm cri}} \kappa (r, \theta, z)\, \rho (r, \theta, z)\,dz}{\sigma _{\rm s}(r, \theta)}}. \end{equation} \tag{ 15 }$$

Zoom In Zoom Out Reset image size

Then, using Equations (12), (14), and (15), we derived the azimuthal-averaged Q parameter as

$$\begin{equation} Q (r) = {\displaystyle \frac{dr \int _0^{2 \pi } Q(r, \theta)\, r_c\, d\theta }{dS}}, \end{equation} \tag{ 16 }$$

where

$$\begin{equation} Q(r, \theta) = {\displaystyle \frac{c_{\rm s,s}(r, \theta)\, \kappa _s (r, \theta)}{\pi G\, \sigma _r (r, \theta)}}. \end{equation} \tag{ 17 }$$

The azimuthal-averaged plasma β_p(≡8πc²_sρ/B²) is also derived in a similar way to the Q parameter as

$$\begin{equation} \beta _{\rm p} (r) = {\displaystyle \frac{dr \int _0^{2 \pi } \beta _{\rm p, s} (r, \theta)\, r_{\rm c}\, d\theta }{dS}}, \end{equation} \tag{ 18 }$$

where

$$\begin{equation} \beta _{\rm p,s}(r, \theta) = {\displaystyle \frac{\int ^{z < z_{\rm cri}} \beta _{\rm p} (r, \theta, z)\, \rho (r, \theta, z)\,dz}{\sigma _{\rm s}(r, \theta)}}. \end{equation} \tag{ 19 }$$

In addition, the ratio of the azimuthal-averaged v_ϕ to the Kepler velocity on the equatorial plane is plotted in Figure 10(d). The black lines in Figure 10 are the radial quantities at the epoch just before the first fragmentation (t_c = 445 yr) for model A03, while the blue lines are those at the epoch just before the second fragmentation after the first fragments fall onto the protostar (t_c = 3947 yr) for model A03. The red-dotted lines are quantities at almost the same epoch (t_c = 3781 yr) as the blue lines but for model A05.

The black lines in Figures 10(a) and (b) show that, for r < 10 AU, the surface density has a maximum and the Q parameter becomes Q < 1. Just after this epoch, fragmentation occurs at r = 6.8 AU where Q < 1. The black lines in Figures 10(c) and (d) show that the plasma β_p suddenly rises at r ∼ 10 AU, inside which the flow moves azimuthally with close to Keplerian velocity. Since the gas (surface) density is sufficiently high in the region of r < 10 AU, the Ohmic dissipation becomes effective and the magnetic field considerably dissipates. In general, fragmentation hardly occurs when the neutral gas is strongly coupled with the magnetic field (or ions), because the magnetic field suppresses fragmentation (Hennebelle & Teyssier 2008; Machida et al. 2008a). Instead, when Q < 1, fragmentation easily occurs in the magnetically inactive region where the Ohmic dissipation weakens the magnetic field. In summary, at this epoch, fragmentation can occur in the region of r < 10 AU, because Q < 1 is realized owing to the higher surface density and the very weak magnetic field (or higher β_p) does not suppress fragmentation.

In the region of r < 10 AU, the surface density at t_c = 3947 yr (blue line) is lower than that at t_c = 445 yr (black line). This is because, in the inner disk region, the angular momentum is transferred by the orbital motion of the first-generation planets and the gas falls onto the protostar. In addition, planets contribute to the gap formation around the central protostar as seen in Vorobyov & Basu (2010). Thus, at t_c = 3947 yr, the surface density near the protostar becomes lower and the Q parameter exceeds Q>1 for r ≲ 30 AU (blue line of Figure 10(b)). Since the lower (surface) density makes the Ohmic dissipation ineffective, a lower plasma β_p (β_p ∼ 2–4) is realized in this region (r ≲ 40 AU). Therefore, the inner disk region becomes magnetically active and the angular momentum is also transferred by magnetic braking. Figure 10(d) shows that the ratio of the azimuthal velocity to the Kepler velocity at t_c = 3947 yr is lower than that at t_c = 445 yr. The lower surface density and higher value of plasma β_p are not adequate to induce fragmentation. On the other hand, the surface density and plasma β_p at t_c = 3947 yr have a local peak at r ∼ 30–40 AU (blue lines of Figures 10(a) and (c)). This is because, in this region, the radially inward motion of the flow is suppressed by the first-generation planet and the gas is accumulated. In such a high-density gas region, Ohmic dissipation becomes effective owing to the higher surface density. Thus magnetic braking becomes ineffective and azimuthal velocity approaches the Kepler speed (blue line of Figure 10(d)). Therefore, the centrifugal barrier further amplifies the surface density. Finally, the Q parameter becomes Q < 1 in the range of 20 AU ≲ r ≲ 60 AU (Figure 10(b)), and fragmentation occurs to form second-generation planets at r ≃ 46 AU as seen in Figure 1(g).

In Figures 10(a) and (b), model A05 (red dotted line) has a relatively lower surface density because this model has a lower mass accretion rate onto the circumstellar disk. Therefore, the Q parameter never reaches Q ≪ 1, and fragmentation does not occur until the end of the calculation for model A05. In addition, the lower surface density increases the magnetically active area (Figure 10(c)) and a stronger outflow appears as seen in Figure 6. Thus, the excess angular momentum is mainly transferred by the protostellar outflow and magnetic braking for model A05, while the excess angular momentum contributes to fragmentation for model A03.

In summary, the initial stability (or mass accretion rate) is related not only to the surface density (or mass of the circumstellar disk) but also to the magnetic activity in the circumstellar disk. The higher accretion rate has a higher surface density (i.e., lower Q) and has a large magnetically inactive area where the angular momentum is not transferred by the magnetic effects. Both effects promote fragmentation and planet formation in the circumstellar disk. Instead, since the lower accretion rate has a lower surface density and large magnetically active area, fragmentation and planet formation tend to be suppressed.

The condition of the circumstellar disk also affects protostellar outflow. Figure 11 shows the mass outflowing from the protostellar system (i.e., protostar + circumstellar disk system) for both models. To measure the outflowing mass, we integrated the gas having v_r>c_s,0 inside the initial radius of the cloud core (r < R_c), where v_r is the radial velocity. The outflowing mass for model A05 gradually increases for t_c ≲ 6000 yr and decreases for t_c ≳ 6000 yr. Since the free-fall timescale of the initial cloud core is ∼1.8 × 10⁴ yr, it is expected that the outflow weakens after passing the peak and disappears on a free-fall timescale. In addition, the outflow extends far beyond the initial radius of the cloud core for t_c ≳ 5000 yr. Thus, we underestimated the outflowing mass for t_c ≳ 5000 yr because we integrated the outflowing mass only inside r < R_c. To estimate conclusively the outflowing mass and total ejected mass from the initial cloud core, more long-term calculations are necessary. However, since such calculations are beyond the scope of this study, we no longer comment on them. At its peak, the outflowing mass is about 0.05 M_☉ which is comparable to the protostellar mass M_ps ≃ 0.1 M_☉ (see Figure 7) at the same epoch. Thus, a large fraction of the infalling gas is blown away from the circumstellar disk by the protostellar outflow for model A05.

Zoom In Zoom Out Reset image size

Unlike model A05, in Figure 11, the outflow mass for A03 shows high time variability, indicating that outflow is unsteadily driven by the circumstellar disk. In summary, Figure 11 indicates that steady outflow is driven by the less massive circumstellar disk without fragments (model A05), whereas intermittent outflow is driven by the more massive circumstellar disk with fragments (model A03). Note that Figure 11 shows only the outflowing mass having the high-velocity component of v>c_s,0, and ignores the low-velocity component of the outflow.

Figure 12 shows the ratio of the outflowing to infalling mass. To investigate the outflow efficiency, we calculated the outflowing/infalling mass rate on an l = 9 grid surface, which covers the whole circumstellar disk with a grid size of ∼600 AU. According to Tomisaka (2002), we estimated the outflow $\dot{M}_{\rm out}$ and inflow $\dot{M}_{\rm in}$ masses as

$$\begin{equation} \dot{M}_{\rm out} = \int _{\rm boundary\, of \, l=9} \rho \, {\rm max}[\mbox{\boldmath $v$} \cdot \mbox{\boldmath $n$},0]\, {\it dS} \end{equation} \tag{ 20 }$$

and

$$\begin{equation} \dot{M}_{\rm in} = \int _{\rm boundary\, of \, l=9} \rho \, {\rm max}[\mbox{\boldmath $v$} \cdot \mbox{\boldmath $-n$},0]\, {\it dS}, \end{equation} \tag{ 21 }$$

where $\mbox{\boldmath $n$}$ is the unit vector outwardly normal to the surface of the l = 9 grid. In Figure 12, model A05 keeps an almost constant ratio, $\dot{M}_{\rm out}/\dot{M}_{\rm in}\sim 0.3$, until the end of the calculation. Thus, about 30% of the infalling mass is ejected by the outflow. In this model, the outflow is driven by the outer disk region where the axisymmetric structure is kept, while no outflow appears in the inner disk region where a considerably weak magnetic field is realized by Ohmic dissipation and a non-axisymmetric structure appears, as seen in Figure 6. In the outer disk region, the disk mass is blown away by the outflow, and effective magnetic braking moves the gas inward. Therefore, such a region does not maintain sufficient mass to induce gravitational instability. As a result, the steady disk region drives the steady flow for model A05.

Zoom In Zoom Out Reset image size

On the other hand, the outer disk region for model A03 does not remain stable, because the increase rate of the disk mass (or the mass accretion rate onto the disk) for model A03 is much larger than that for model A05. Although magnetic braking and mass ejection by the outflow are effective in the outer disk region for model A03 also, they cannot transfer sufficient mass to suppress gravitational instability and subsequent fragmentation owing to the high mass accretion onto the disk. Thus, multiple planets appear in the disk. Oscillation of the infall–outflow mass ratio for model A03 is shown in Figure 12. Comparison of Figure 4 with Figure 12 indicates that this oscillation is related to the orbital motion of the planet. It is considered that the magnetic field lines and surface density of the circumstellar disk are disturbed by the planet's orbital motion, which causes a highly time-variable outflow, as shown in Figure 12. The time-variable protostellar outflow may be a clue for detecting a planet in a very early phase of star formation.

So far, it is considered that planets form in relatively quiet circumstellar disks after the gas accretion onto the circumstellar disk has almost ended. In this study, however, we showed the possibility of the formation of planets or substellar-mass objects by gravitational instability in a vigorous circumstellar disk during the main accretion phase, in which the circumstellar disk continues to increase its mass by gas accretion from the infalling envelope. When planets or substellar-mass objects appear in the circumstellar disk, they can affect the rate of accretion onto the protostar. As described in Section 3.5, we may verify their existence by observing the time variability of the rate of accretion onto the central protostar. As seen in Figure 9, the protostar shows a relatively low time variability of the accretion rate when neither planets nor substellar-mass objects exist in the circumstellar disk, while it shows a high time variability of the accretion rate when fragmentation occurs and planets or substellar-mass objects appear. Figure 13 shows the time variation of the accretion luminosity of the protostar, which is defined as

$$\begin{equation} L_{\rm ps} = G\, {\displaystyle \frac{M_{\rm ps}\, \dot{M}_{\rm ps}}{r_{\rm ps}}}, \end{equation} \tag{ 22 }$$

where we adopted the protostellar mass M_ps and mass accretion rate $\dot{M}_{\rm ps}$ from values derived from the calculations at each time and used a constant protostellar radius r_ps = 2 R_☉ for simplicity. The figure indicates that the protostar for model A05 has an accretion luminosity in the range of $0.01 \lesssim L_{\rm ps}/\thinspace L_\odot \lesssim 1$ for ∼10⁴ yr after protostar formation. The accretion luminosity for this model gradually increases for a longer timescale of ∼10⁴ yr, while it oscillates with a period of ∼100 yr, which corresponds to the Kepler timescale in the region of a developing spiral structure. The Kepler timescale τ_k can be expressed as

$$\begin{equation} \tau _{\rm k} = 2\pi \left({{\displaystyle \frac{r^3}{{\it G M}_{\rm p}}} } \right)^{1/2} \simeq 100 \left({\displaystyle \frac{r}{10\,{\rm AU}}} \right)^{3/2} \left({\displaystyle \frac{0.1\, M_\odot }{M_{\rm ps}}} \right)^{-1/2}\, {\rm yr}. \end{equation} \tag{ 23 }$$

Zoom In Zoom Out Reset image size

As shown in Figure 6, the non-axisymmetric structure appears at r ≲ 10–40 AU, which corresponds to the magnetic inactive zone (Machida et al. 2007). Thus, it is natural that the protostar shows a time-variable accretion luminosity with a period of ∼100 yr, because the angular momentum transfer that is closely related to the rate of mass accretion onto the protostar is caused by the spiral structure in the circumstellar disk. The variation in the amplitude of the accretion luminosity for a short duration (∼100 yr) is about one order of magnitude.

On the other hand, the accretion luminosity for model A03 is qualitatively and quantitatively different from that for model A05. The protostar for model A03 repeatedly shows a strong brightening. In this model, the accretion luminosity suddenly increases to reach ${\sim} 10\hbox{--}10^4\thinspace L_\odot$ when planets or substellar-mass objects approach the protostar or fall onto the protostar, while the luminosity is lower, ≲0.1–0.01 L_☉, when planets or substellar-mass objects orbit far from the protostar. Thus, the strong brightness may be a clue for detecting planets or sub-stellar-mass objects.

In addition, when the planets or substellar-mass objects orbit in the inner disk region of ≲10–20 AU, their orbital motion is linked to the accretion rate (or accretion luminosity). The accretion rate and orbital radius for t_c = 1000–2800 yr for model A03 are plotted in Figure 14. In this time period, the first-generation planet orbits around the protostar for ∼10 times. The figure indicates that the orbital motion of planets or substellar-mass objects clearly synchronizes with the mass accretion rate. The orbiting objects disturb the circumstellar disk and promote mass accretion onto the protostar. Since the circumstellar disk has a higher density and shorter timescale in the proximity of the protostar, a high mass accretion rate is realized when the planets or substellar-mass objects approach the protostar. As a result, the mass accretion rate is closely related to the orbital motion of planets and substellar-mass objects. We may estimate and verify the orbital motion of planets in the circumstellar disk by observing the time-variable mass accretion. Note that we ignored the effect of the infalling envelope when the accretion luminosity is estimated in Figure 13. Note also that, to estimate quantitatively the accretion luminosity, we have to consider the infalling envelope that weakens the luminosity from the protostar.

Zoom In Zoom Out Reset image size

In this study, we showed that fragmentation and subsequent formation of planet and sub-stellar mass object can occur in a massive circumstellar disk, while no fragmentation occurs in a relatively less massive disk. The mass of the circumstellar disk or the mass increase rate of the circumstellar disk is related to the initial stability of cloud core. We parameterized the ratio of thermal energy to gravitational energy, α₀, which is related to the mass accretion rate. In general, the mass accretion rate is described as

$$\begin{equation} \dot{M} = {\displaystyle \frac{M_{\rm J}}{t_{\rm ff}}} \propto {\displaystyle \frac{c_s^3}{G}}, \end{equation} \tag{ 24 }$$

where M_J is the Jeans mass and t_ff is the free-fall timescale. In this study, we made an equilibrium sphere (i.e., critical BE sphere) and increased the density (or mass) by a factor of f (=0.84/α₀, see Equation (6)) to induce the collapse of the cloud core. Thus, the mass accretion rate in our setting can be described as

$$\begin{equation} \dot{M} = {\displaystyle \frac{M_{\rm cl}}{t_{\rm cl,ff}}} = \alpha _0^{-3/2}\, {\displaystyle \frac{c_s^3}{G}}, \end{equation} \tag{ 25 }$$

where the initial mass of the cloud core M_cl is proportional to the density enhancement factor f and Jeans mass M_J as M_cl ∝ f M_J ∝ α⁻¹ M_J,and the free-fall timescale of the cloud core t_ff,cl is proportional to t_ff,cl ∝ f^−1/2 t_ff,0 ∝ α^1/2 t_ff,0,where t_ff,0 is the free-fall timescale of a critical BE sphere. Thus, a more unstable cloud core that has a smaller α₀ has a larger mass accretion rate and tends to show fragmentation due to gravitational instability. In reality, the mass increase rate for the model A03 is about three times higher than that for the model A05. Thus, the fragmentation and subsequent planets or substellar-mass objects may appear in relatively unstable cloud cores.

In general, however, the initial clouds are characterized by three parameters, α₀, β₀, and γ₀, as described in Section 2.2. Although we fixed the angular velocity (β₀) and magnetic field strength (γ₀) of the initial cloud core in this study, fragmentation may occur in more stable cloud cores (i.e., cloud cores with large α₀) with different β₀ and γ₀. As shown in Hennebelle & Teyssier (2008), in the collapsing cloud core, the magnetic field suppresses fragmentation while the rotation promotes it. Thus, fragmentation tends to occur in the cloud core with larger β₀ and smaller γ₀ (Machida et al. 2008b). Moreover, even in a cloud core with a very weak magnetic field, fragmentation in the circumstellar disk may be suppressed by the global spiral structure that effectively transports the angular momentum and lowers the surface density (i.e., increases the Toomre Q-parameter). On the other hand, when the circumstellar disk is strongly magnetized, the angular momentum is transferred by magnetic effects and the spiral structure hardly develops. Then, the mass accumulates in the inner disk region where Ohmic dissipation significantly weakens the magnetic field, and fragmentation can occur because no global spiral structure develops. Thus, magnetic field and rotation contribute to both promotion and suppression of fragmentation through a complicated process. In this study, we showed the possibility of fragmentation in the main accretion phase. However, to quantitatively determine fragmentation and the subsequent formation of planets and substellar-mass objects in the circumstellar disk, we should investigate the evolution of the circumstellar disk from the molecular cloud core in a large parameter space.

In the circumstellar disk, fragmentation occurs and planets or substellar-mass objects appear. In this study, we resolved the fragmentation process with sufficient spatial resolution (Truelove et al. 1997). However, higher spatial resolution and further developed numerical techniques may be necessary to determine the final mass of fragments. At its formation, a fragment's mass almost corresponds to the Jeans mass. In this study, we assumed the barotropic equation of state, in which the minimum Jeans mass is limited to ∼5 M_Jup. Thus, fragmentation does not occur with a fragment mass of M ≲ 5 M_Jup. If the circumstellar disk cools and reaches lower temperatures, more less-massive fragments may appear. The barotropic approximation makes it possible to calculate the very long-term evolution of the circumstellar disk and to investigate the evolution of the circumstellar disk by the end of the main accretion phase from the prestellar core stage, while radiation hydrodynamic calculations are necessary to qualitatively estimate the fragment mass.

In addition, after fragmentation occurs in the circumstellar disk, we need a higher spatial resolution to investigate the acquisition process of the fragments' mass. Recent works about a gas-giant planet formation showed that the Hill radius of the protoplanet should be resolved with sufficient spatial resolution for an accurate estimate of the rate of mass accretion onto the protoplanet (Kley et al. 2001; D'Angelo et al. 2003; Bate et al. 2003; Machida et al. 2010b). Note that in these calculations, the protoplanet mass was fixed, and the evolution of the circumstellar disk was calculated with a high spatial resolution for a shorter duration after a gas giant is formed. The Hill radius is expressed as

$$\begin{equation} r_{\rm H} = \left({\displaystyle \frac{M_{\rm p}}{3\, M_{\rm ps}}} \right)^{1/3} r_{\rm p}, \end{equation} \tag{ 26 }$$

where r_p is the orbital radius of the planet. When the protoplanet orbiting the solar-mass protostar has a mass of M_p = 10⁻³M_ps (10⁻²M_ps) and is located at r = 1 AU (10 AU), the Hill radius is r_H ≃ 0.07 AU (1.5 AU). When the Hill radius is resolved with ≳10 cells, we need a minimum spatial resolution (i.e., the size of each cell) of ∼0.007 AU (∼0.15 AU). However, calculation with such a high resolution needs a lot of computing power, and thus we only calculated the evolution of the circumstellar disk for a very short duration. In reality, the evolution of the circumstellar disk was calculated for only ∼1000 yr at best in previous studies, while we calculated the evolution of the circumstellar disk from the prestellar core stage until ∼10⁴ yr after protostar formation.

The final mass of the fragments is important to determine whether they are planets or binary companions. To understand both formation and evolution of (gas-giant) planets, we need to calculate the evolution of the circumstellar disk for a longer duration with a higher spatial resolution. However, such calculations are impossible even using present supercomputers because a longer temporal calculation conflicts with higher spatial resolution. Thus, in this study, we covered the whole prestellar cloud core and investigated the evolution of the circumstellar disk with moderately coarser spatial resolution Δ ≃ 0.6 AU, and showed that fragmentation occurs in the circumstellar disk and a planet or substellar-mass object forms. However, since we need a higher spatial resolution to determine the final mass of the fragment, we did not investigate this point in greater depth. In this study, we showed the possibility of the formation of planets or substellar-mass objects in a magnetized disk for the main accretion phase, and the precise mass of such objects should be investigated in future studies.

As described in Sections 3.6 and 3.7, the recurrent planet formation and intermittent outflow shown in this paper are closely related to the effect of the magnetic field and its dissipation in the circumstellar disk. Figure 10 shows that the magnetic field is well coupled with the neutral gas in the outer disk region (r_c>10–100 AU) where a relatively high ionization rate is realized owing to a relatively low gas density. Thus, the angular momentum in such region is effectively transferred by the magnetic braking and protostellar outflow, and the gas can flow into the inner disk region. On the other hand, the magnetic field is not well-coupled with the neutral gas in the inner disk region (r_c < 10–100 AU) where the gas density is high and a considerably low ionization rate is realized. Thus, the angular momentum transfer due to the magnetic effects is not so effective, and the gas flowing from the outer disk region accumulates in the inner disk region around the boundary between the magnetically active and inactive region at r_c ∼ 5–50 AU (see Figure 10). In that region, the surface density continues to increase and fragmentation tends to occur. In summary, magnetic dissipation is essential in inducing recurrent fragmentation (or planet formation) in the circumstellar disk.

In this study, we have considered only Ohmic dissipation as the dissipation mechanism of the magnetic field. This is because Ohmic dissipation is expected to be the most effective mechanism for magnetic field dissipation in a high density region where most of the magnetic flux leakage occurs. However, ambipolar diffusion may also contribute to the magnetic flux leakage in some intermediate density region. Recent studies have shown that Ohmic dissipation is more effective than ambipolar diffusion in the range of n ≳ 10¹²  cm⁻³, while ambipolar diffusion dominates in the range of 10¹¹  cm⁻³ ≲ n ≲ 10¹²  cm⁻³ (Tassis & Mouschovias 2007a, 2007b, 2007c; Kunz & Mouschovias 2010). In the range of n ≲ 10¹¹  cm⁻³ the magnetic field is well-coupled with the neutral gas in a magnetically super critical cloud (Nakano et al. 2002; Kunz & Mouschovias 2010). Thus, it might be important to study how the inclusion of ambipolar diffusion would change the formation and evolution of the circumstellar disk, through the possible contribution in the limited range of density.

With strong coupling approximation, the induction equation including the ambipolar diffusion can be written as

$$\begin{eqnarray} {\displaystyle \frac{\partial \mbox{\boldmath $B$}}{\partial t}} &=& \nabla \times (\mbox{\boldmath $v$} \times \mbox{\boldmath $B$}) + \eta _{\rm OD} \nabla ^2 \mbox{\boldmath $B$} + \nabla\nonumber\\ && \times \left[ {\displaystyle \frac{\mbox{\boldmath $B$}}{4\pi \gamma \rho _n \rho _i}} \times \left[ \mbox{\boldmath $B$} \times (\mbox{\boldmath $\nabla $} \times \mbox{\boldmath $B$}) \right] \right], \end{eqnarray} \tag{ 27 }$$

where ρ_n, ρ_i, and γ are the densities of neutral gas and charged particles, and a drag coefficient between charged particles and neutrals, respectively (for details, see, e.g., Shu 1992). The third term of the right-hand side of Equation (27) can be regarded as a nonlinear diffusion term by the ambipolar diffusion. When we adopt an approximate relation between the neutral and charged particle densities, ρ_i = Cρ^1/2_n, where C is a constant, the effective diffusion coefficient of the ambipolar diffusion (η_ad) can be roughly estimated as

$$\begin{equation} \eta _{\rm ad} = {\displaystyle \frac{B^2}{4\pi \gamma C \rho ^{3/2}}}. \end{equation} \tag{ 28 }$$

In a realistic setting, Kunz & Mouschovias (2010) showed that the diffusion coefficient of the ambipolar diffusion (η_ad) is 10–100 times larger than that of the Ohmic dissipation (η_OD) in the range of 10¹¹  cm⁻³ ≲ n ≲ 10¹²  cm⁻³. Strictly speaking, to include the effect of ambipolar diffusion, we should solve Equation (27) together with the equations determining charged particle density. Note that the realistic evolution of charged particle density sensitively depends on the total surface area of dust grains, since the recombination of charged particles on the grain surface is very efficient. However, there is considerable uncertainty in the evolution of dust grain properties, such as the size distribution, in the protostellar collapse process. Thus, our limited knowledge of grain properties make our accurate treatment of magnetic field dissipation impractical. This is also the reason why we adopted a parametric description of the Ohmic diffusion coefficient in our study of magnetic field dissipation (Machida et al. 2007). In this paper, we take the same approach for the effect of ambipolar diffusion. To roughly estimate the effect of ambipolar diffusion, we calculate the evolution of the circumstellar disk with an artificially larger resistivity. In Equation (3), we replace η_OD by c_η η_OD as follows:

$$\begin{equation} {\displaystyle \frac{\partial \mbox{\boldmath $B$}}{\partial t}} = \nabla \times (\mbox{\boldmath $v$} \times \mbox{\boldmath $B$}) + c_\eta \, \eta _{\rm OD} \nabla ^2 \mbox{\boldmath $B$}. \end{equation} \tag{ 29 }$$

We adopt c_η = 0 (ideal MHD model), 1 (fiducial model), 10, 100.

As a first test, we constructed a non-rotating BE sphere with a central density of n ≃ 10³  cm⁻³ immersed in the uniform magnetic field of B_z,0 = 5 × 10⁻⁵ G. Then, we calculated the collapse of this cloud with different c_η. Figure 15 shows the evolution of the magnetic field as a function of the central number density. The magnetic field in each model tracks the same path in the range of n ≲ 10⁹  cm⁻³, while the magnetic field for the model with c_η = 100 begins to depart from that for the ideal MHD model (c_η = 0) at n ≃ 10⁹  cm⁻³. The figure indicates that the non-ideal MHD effect becomes important at the lower density with the larger c_η. When the central number density reaches n = 10¹⁴  cm⁻³, the magnetic field strengths for the model with c_η = 0 and 100 are B_z = 7 G and 0.03 G, respectively. Thus, the ideal MHD model has about a 100 times stronger magnetic field than the model with c_η = 100 at this epoch. At the same epoch, models with c_η = 1 and 10 have B_z = 0.01 G and 0.06 G of magnetic field, respectively. Thus, the difference between the magnetic field strength between c_η = 1 and 100 (10) is only a factor of three (2).

Zoom In Zoom Out Reset image size

Next, we calculated the evolution of the circumstellar disk with the initial condition the same as model A03 but with different resistivities c_η = 0, 10, and 100. The density distribution on the equatorial plane after the circumstellar disk formation for models with c_η = 0 and c_η = 100 are plotted in Figure 16. As seen in this figure, fragmentation does not occur in the ideal MHD model (c_η = 0) until the end of the calculation, while, for the model with c_η = 100, fragmentation occurs at t_c ≃ 2550 yr after the protostar formation and two clumps appear in the circumstellar disk. For the model with c_η = 10, the circumstellar disk also shows fragmentation and subsequent planet formation. Note that we have terminated the calculation with larger resistivity earlier, since it takes formidable time to follow the evolution with larger resistivity with our time-explicit integration scheme for diffusion term (Machida et al. 2007). Thus, we cannot directly compare the results among models with different resistivities. However, Figure 16 indicates that fragmentation occurs even when the magnetic field dissipates in a relatively low-density region. Therefore, we expect that recurrent planet formation and intermittent outflow can occur even with a more realistic treatment of the ambipolar diffusion.

Zoom In Zoom Out Reset image size

Only a few studies focused on the formation and evolution of the circumstellar disk from the prestellar core stage. Walch et al. (2009a) and Vorobyov & Basu (2010) investigated the circumstellar disk in unmagnetized clouds, and found that massive circumstellar disks formed at the early main accretion phase tend to show fragmentation. Machida et al. (2010a) and Bate (2010) pointed out that the disk forms before the protostar formation, and the circumstellar disk is inevitably more massive than the protostar in the early accretion phase. This is because the circumstellar disk originated from the first adiabatic core formed in the gas-collapsing phase before the protostar formation (Bate 1998; Machida et al. 2010a). Such a massive circumstellar disk is prone to show fragmentation. These studies indicate the possibility of the formation of a substellar- or planet-mass object in the unmagnetized circumstellar disk by gravitational instability.

In the magnetized cloud core, only Vorobyov & Basu (2005, 2006) investigated the formation and evolution of the circumstellar disk using two-dimensional simulations. Although they adopted the magnetic evolution in an approximate way instead of solving the induction equation, they presented several interesting results. They pointed out that fragmentation and planet formation is possible even in the magnetized circumstellar disk. They also showed that gravitational instability or fragments temporarily amplify the mass accretion rate onto the protostar that may account for the FU Orionis outbursts. There are slight quantitative differences between Vorobyov & Basu (2005, 2006) and our results, because we calculated the circumstellar disk with a realistic process of the magnetic dissipation in three dimensions. However, our results are qualitatively the same as their results. This is because the circumstellar disk always experiences a gravitationally unstable phase after protostar formation. Inutsuka et al. (2010) showed that, even in the magnetized cloud core, the circumstellar disk is more massive than the protostar and becomes gravitationally unstable. Such a disk tends to show fragmentation in the early accretion phase. Even when fragmentation does not occur in the early accretion phase, such a massive disk is expected to show fragmentation in the later accretion or subsequent phase. Thus, we expect that the planet formation due to the gravitational instability is common in the star formation process.