In our present view of low-mass star formation, a protostar in the early stages of mass assembly (in the so-called class 0 and class I phases) is surrounded by a protostellar disk which is in turn deeply embedded in an infalling envelope left over from the collapse of a rotating prestellar cloud core. The observed low luminosity of these protostars implies a low mass accretion rate, and hence a long lifetime in order to achieve typical final stellar masses; however, this is inconsistent with the number of known class 0 and class I objects (Kenyon et. al. 1990). A possible explanation is that the protostellar accretion proceeds in two co-existing phases (Kenyon & Hartmann 1995). Accretion from the envelope onto the protostellar disk takes place in a fairly uniform (though generally declining in time) manner, whereas accretion from the protostellar disk onto the central protostar occurs primarily in short (and infrequently observed) but powerful episodes during which 0.01−0.1 M_⊙ can be accreted. These episodes of vigorous accretion

≥ 10⁻⁴ M_⊙ yr⁻¹ manifest themselves as FU Orionis variables (FU Ori). Between these accretion bursts, a typical class 0/class I protostar is characterized by a low accretion rate

∼ 10⁻⁷ M_⊙ yr⁻¹.

The nature of FU Ori disk accretion bursts has been widely debated. For instance, close encounters in binary systems may cause a strong perturbation in protostellar disks and drive high accretion rates during a relatively short period Bonnell. This mechanism requires rather eccentric orbits of binary systems and obviously fails to explain FU Ori outbursts in isolated protostars. An alternative idea is that the thermal instability (namely, the steep dependence of the disk opacity on temperature between ∼ 3000 K and ∼ 10⁴ K) of the optically thick innermost regions of circumstellar disks triggers FU Ori eruptions (Lin & Papaloizou 1985; Clarke et al. 1990; Bell & Lin 1994). Most thermal instability models exploit the α-prescription of Shakura, who suggested that the disk effective viscosity is proportional to its temperature; an increase in disk temperature causes a higher rate of mass accretion due to an elevated viscous mass transfer and vice versa. Unfortunately, many aspects of this mechanism of FU Ori outbursts are completely dependent upon the unknown value of disk effective viscosity, which determines the timescale for outbursts.

It is known from theoretical and numerical studies that protostellar disks may be subject to the development of global spiral instabilities. For instance, Laughlin & Bodenheimer (1994) have studied the nonaxisymmetric evolution of protostellar disks and found that they are susceptible to a series of spiral instabilities, with the fastest growing modes being the one-armed (m=1) and two-armed (m=2) patterns. Recent numerical simulations of star cluster formation (Bate et al. 2003) confirmed that the protostellar disks formed around protostars were gravitationally unstable and prone to the development of spiral density waves. Simulations of marginally unstable protoplanetary disks show the formation of flocculent and clumpy spiral structure, suggesting a possible transient rise in the mass accretion rate associated with clump infall (Boss 2003). Mejia et al. (2005) have also reported three-dimensional simulations which demonstrate a single FU-Ori-like outburst associated with the growth of spiral structure in an isolated protoplanetary disk. Recent observations do reveal a flocculent spiral structure in the protostellar disk of AB Aurigae (Fukagawa et al. 2004; Corder et al. 2004).

In this Letter, we present the first model of cloud core collapse which self-consistently generates multiple accretion bursts that can be identified with FU Ori eruptions. The protostellar disk in our model is formed as a result of the collapse and is not isolated from the parent core envelope. The details of the numerical model are described in § 2. The results of simulations are presented in § 3. Our main conclusions are summarized in § 4.

2. Model description

We model the collapse of a rotating cloud core which is threaded by a frozen-in magnetic field with spatially uniform mass-to-flux ratio. The magnetic field effect is comparable to but weaker than gravity, so that the core is magnetically supercritical. A magnetically-diluted gravitational collapse ensues. Our initial core model is a good approximation to the supercritical cores that result from the fragmentation of magnetized clouds that are initially either critical/subcritical or supercritical (e.g., Basu 1997; Basu & Ciolek 2004).

We follow the collapse through the prestellar collapse phase and into the accretion phase which sees the development of a protostar and a protostellar disk. The thin-disk approximation is used in the form appropriate for a supercritical core with mass-to-flux ratio μ that is spatially uniform (Basu 1997; Nakamura & Hanawa 1997; Shu & Li 1997). The magnetic field pressure enhances the gas pressure and the magnetic tension dilutes the effect of gravity. The disk is symmetric about the midplane, with a vertical magnetic field component B_z inside the disk and both vertical and tangential components outside it. The material external to the disk is current-free (j = 0). We characterize the field strength by the parameter α ≡ μ⁻¹ = B_z/ (2 π √GΣ), where Σ is the gas surface density; note α < 1 for a supercritical core.

The magnetohydrodynamic equations we use are written in the thin-disk approximation as

We assume a two-component equation of state P=c_s²Σ + c_s² Σ_cr (Σ/Σ_cr)^γ, where c_s is the isothermal speed of sound and Σ_cr is the critical gas surface density above which the disk becomes optically thick. This equation of state allows for a smooth transition between the isothermal and non-isothermal regimes during the collapse. We use a canonical value of the critical gas volume density n_cr=10¹¹ cm⁻³ (Larson 2003), which is equivalent to Σ_cr=36.2 g cm⁻² for the gas disk in vertical hydrostatic equilibrium. In equation (2), we assume a cloud scale height Z=c_s²/(π G Σ) for Σ≤ Σ_cr and Z=c_s²/(π GΣ_cr) for Σ > Σ_cr. We adopt the ratio of specific heats γ = 7/5 for the optically thick regime, appropriate for an adiabatic diatomic gas.

Equations (1) and (2) are numerically solved in polar coordinates (r,φ) using the method of finite differences with a time-explicit, operator-split solution procedure similar to that described by Stone and Norman in their ZEUS-2D code (Stone & Norman 1992). The details of the code and results of the tests will be given in a follow-up paper. The numerical grid has 256×256 points, which are logarithmically spaced in r-direction allowing for a good resolution in the inner cloud regions. The innermost grid point is located at 10 AU and the size of the first adjacent cell is 0.3 AU. The Truelove criterion (Truelove et al. 1998) is preserved throughout the simulations – the size of grid cells is smaller than the Jeans length. We introduce a “sink cell” at r < 10 AU, which represents the central protostar plus some circumstellar disk material, and impose a free inflow inner boundary condition. We assume that the matter is cycled through the circumstellar disk and onto the protostar rapidly enough so that the mass infall through the sink cell is at least proportional to the mass accretion rate onto the protostar. We impose the outer boundary condition such that the gravitationally bound cloud has a constant mass and volume. The gravitational potential of the cloud is evaluated numerically using the fast Fourier transform (Binney & Tremaine 1987).

3. Results

The cloud evolution is characterized by a slow initial gravitational contraction and then a very rapid runaway collapse until the formation of the central protostar. The black line in Figure 1a shows two distinct phases in the temporal evolution of the mass accretion rate onto the protostar. The early behavior of is qualitatively similar to that obtained in spherical collapse simulations (see e.g., Vorobyov & Basu 2005). Accretion shows a very rapid increase to a maximum value of =1.0 10⁻⁴ M_⊙ yr⁻¹ at t = 0 yr, when the central protostar forms. Subsequently, there is a slow decline in , when the gas is accreted directly onto the protostar from the inner envelope, which has a relatively low specific angular momentum. The second phase starts at t = 7700 yr when the protostellar disk forms around the protostar due to the infall of matter from the outer envelope, which has a higher specific angular momentum. The mass of the protostar at this stage is approximately M_s = 0.4 M_⊙. The matter in the disk moves on nearly circular orbits and most of it is far from the protostar. Hence, the accretion rate onto the protostar abruptly drops down to a negligible value. A continuous infall of matter from the envelope makes the protostellar disk unstable to the development of spiral structure shown in Fig. 2 and induces the formation of dense protostellar/protoplanetary clumps within the arms. The even (m = 2,4) modes are dominant in this simulation, but the odd (m = 1,3) modes are also excited in some simulations. Spiral arms transport angular momentum outward and mass inward (Lynden-Bell & Kalnajs 1972). The ensuing redistribution of mass and angular momentum creates a strong centrifugal disbalance in the protostellar disk and triggers bursts of mass accretion when dense protostellar/protoplanetary clumps in the inner disk are driven into the protostar (see below for an animation of the disk evolution, or go here to download the animation).

Animation of the disk evolution.

During this process, the mass in the protostellar disk remains somewhat less than the mass of the protostar. The episodes of clump infall can manifest themselves as very short (≤ 100 yr) but vigorous (

= [1−10] × 10⁻⁴ M_⊙ yr⁻¹) accretion bursts as is clearly seen in Fig. 1. During the accretion bursts, 0.01−0.05 M_⊙ of gas is accreted and the accretion luminosity may by grow many orders of magnitude. The duration of the intervening quiescent accretion phase with

= (1−10) × 10⁻⁷ M_⊙ yr⁻¹ is usually (1−3) × 10³ yr. The frequency of bursts decreases with time and the number of bursts may amount to 15-30. The brown line in Fig. 1b shows a luminosity L_bol = G M_s

/R_c (it is assumed to derive entirely from the disk accretion), where R_c = 4 R_⊙ is the radius of the protostar (Masunaga & Inutsuka 2000). Our L_bol is an upper limit to the expected observable bolometric luminosity. It is evident that the episodes of clump infall lead to a dramatic increase in L_bol (by a factor of up to 2000 as compared to a quiescent period) and the protostar may reveal itself as an FU Orionis variable.

Our calculations of the approximate Toomre parameter Q = c_s Ω/(πGΣ) (Toomre 1981) and the total gravitational torque per unit mass Γ, which is the sum of the individual torques per unit mass (|∂Φ/∂φ|) on all computational cells, support the above scenario of episodic accretion. The quantity c_s² ≡ dP / dΣ is the squared effective sound speed. The Q parameter may serve as an approximate stability criterion – gas disks are gravitationally unstable to local nonaxisymmetric perturbations if Q ≤ 1.5−1.7 (Nelson et al. 1998; Boss 1998), while Γ may roughly express the efficiency of angular momentum and mass redistribution by spiral inhomogeneities in the disk (Tomley et al. 1994). The Toomre parameter is calculated by averaging c_s, Ω, and Σ over all computational cells. The red and blue lines in Fig. 1a and 1b show the evolution of Q and the normalized gravitational torque Γ/Γ_max after the protostellar disk formation, respectively. In the early phase of near constant accretion, the matter is directly accreted onto the protostar and Q is much larger than unity. When the protostellar disk forms, its density starts to grow due to accretion. As a consequence, the Toomre parameter gradually decreases below the stability limit Q ≈ 1.5 and reaches a minimum value at the time of the accretion burst. This strongly suggests a causal link between the gravitational instabilities and accretion bursts. The behavior of Γ also shows a direct correlation with accretion bursts. The gravitational torque gradually increases and reaches a maximum value at the time of each accretion burst, indicating the growing efficiency of inward mass transport before the burst. We note that the strength of the torque due to artificial viscosity (used in our Eulerian code to smooth shocks) is at least an order of magnitude smaller than the strength of the gravitational torque associated with spiral instabilities. Thus, the artificial viscosity cannot be responsible for the accretion bursts. After the burst, the mass of the protostellar disk decreases and the growth of spiral instabilities becomes temporarily suppressed, as indicated by high values of Q (> 1.5) and a sharp decrease in Γ. However, the continuous mass infall onto the disk from the envelope quickly destabilizes the disk and the cycle repeats until most of the envelope mass is accreted by the protostar.

Ambipolar diffusion, while not included in our model, is expected to favor the formation of clumps and subsequent burst activity by removing magnetic support. Indeed, a model with no magnetic support (α = 0) shows an increase in the frequency and amplitude of accretion bursts (Vorobyov & Basu, in preparation). Further magnetic effects such as magnetic braking (Krasnopolsky & Königl 2002) and magnetorotational instability (Fromang et al. 2005) may become important in the late accretion phase, but can only be studied in a future three-dimensional model.

4. Conclusions

Rotating protostellar cores show two distinct phases in the temporal evolution of the mass accretion rate

onto the protostar. The early behavior of

is qualitatively similar to that obtained in spherical collapse simulations. Accretion shows a very rapid increase to a maximum, when the central protostar forms, and a subsequent slow decline, when the gas is accreted directly onto the protostar from the inner envelope. The second phase starts when the protostellar disk forms around the protostar due to the infall of matter from the outer envelope with higher specific angular momentum. In this phase,

is characterized by very short (< 100 yr) but vigorous (

= [1−10] × 10⁻⁴ M_⊙ yr⁻¹) accretion bursts, which are intervened with longer periods (~ 10³ yr) of quiescent accretion.

We have demonstrated that the repeating accretion bursts reflect a basic self-regulation mechanism that is inherent to self-gravitating rotating protostellar disks. We emphasize that it is the ongoing infall of matter from the protostellar envelope that continually destabilizes the disk and causes it to periodically dump significant amounts of matter onto the protostar while transferring excess angular momentum to the envelope. The effect of additional support due to a frozen-in supercritical magnetic field is to moderate the burst activity but not suppress it. Recent high resolution Hubble Space Telescope observations of the central regions (≤ a few hundred pc) of Seyfert galaxies (Regan & Mulchaey 1999) also reveal rich spiral structures. We suggest that a self-regulation mechanism similar to that in our model may operate on galactic scales and be responsible for periodic nuclear activity in at least some Seyfert galaxies.

This research was supported by the Natural Sciences and Engineering Research Council of Canada. EIV gratefully acknowledges support from a CITA National Fellowship.

THE ORIGIN OF EPISODIC ACCRETION BURSTS IN THE EARLY STAGES OF STAR FORMATION

E. I. Vorobyov^1,2, Shantanu Basu¹

1. Introduction

2. Model description

3. Results

4. Conclusions

References

THE ORIGIN OF EPISODIC ACCRETION BURSTS IN THE EARLY STAGES OF STAR FORMATION

E. I. Vorobyov1,2, Shantanu Basu1

1. Introduction

2. Model description

3. Results

4. Conclusions

References

E. I. Vorobyov^1,2, Shantanu Basu¹