EDP Sciences logo
Subscriber Authentication Point
Search
Menu
Home All issues Volume 557 (September 2013) A&A, 557 (2013) A35 Full HTML
Free Access
Issue
A&A
Volume 557, September 2013
Article Number A35
Number of page(s) 11
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201321775
Published online 23 August 2013

© ESO, 2013

1. Introduction

Low-mass stars are born as a result of the gravitational collapse of dense and cold cores composed mainly of molecular hydrogen, which in turn are formed in giant molecular clouds from the combination of turbulent, gravitational and magnetic processes (e.g. McKee & Ostriker 2007). Despite many observational and theoretical efforts, how a low-mass star accumulates its final mass remains an open question. In the classic model of inside-out collapse (Shu 1977), the mass accretion rate onto the protostar is proportional to the cube of the sound speed, , implying a rather narrow spread of accretion rates (2 − 5) × 10-6  M⊙ yr-1 for typical conditions in pre-stellar cores. These values of Ṁ are in stark disagreement with those inferred for young protostars, exhibiting a variety of accretion rates from ≲10-7  M⊙ yr-1 to ≳10-5  M⊙ yr-1 (Dunham et al. 2006; Evans et al. 2009). In particular, the fraction of objects with Ṁ < 10-6  M⊙ yr-1 amounts to 50% in Perseus, Serpens, and Ophiuchus star-forming regions (Enoch et al. 2009). FU-Orionis-type objects (FUors) with accretion rates sometimes exceeding 10-4  M⊙ yr-1 and very-low-luminosity objects (e.g. Bourke et al. 2006) with accretion rates most certainly below 10-6  M⊙ yr-1 also do not fit into the classic model.

It is becoming evident that the simplified model of Shu (1977) cannot explain the whole variety of inferred accretion rates and the concept of variable accretion with episodic bursts, based on the original idea of Kenyon et al. (1990), is currently gaining theoretical and observational support. According to this concept, the mass accretion onto the forming protostar is characterized by short (≲100 − 200 yr) bursts approximately with Ṁ ≳ 10-5  M⊙ yr-1 alternated with longer (103 − 104 yr) quiescent periods with Ṁ ≲ 10-6  M⊙ yr-1.

Numerous mechanisms that can produce variable accretion with episodic bursts have been proposed in the past. They include models explaining the origin of FUor accretion bursts by means of viscous-thermal instabilities in the inner disk (Lin & Papaloizou 1986; Bell & Lin 1991), thermal instabilities induced by density perturbations due to, e.g., a massive planet in the disk (Lodato & Clarke 2004), and tidal effects from close encounters in binary systems or stellar clusters (Bonnell & Bastien 1992; Pfalzner et al. 2008). Recently, several promising accretion burst mechanisms have emerged or received further refinement, such as a combination of gravitational instability and the triggering of the magnetorotational instability (Armitage et al. 2001; Zhu et al. 2010), accretion of dense clumps in a gravitationally fragmenting disk (Vorobyov & Basu 2006, 2010), perturbations caused by a planetary or sub-stellar object on an eccentric orbit (Machida et al. 2011a; Vorobyov 2012), or instability when an accretion disk is truncated by the star’s strong magnetic field close to the corotation radius (D’Angelo & Spruit 2010).

Episodic accretion may have important consequences for the evolution of stars and planets. It provides an explanation for the long-standing luminosity problem in young protostars (Dunham & Vorobyov 2012), explains the luminosity spread in young stellar clusters without invoking any significant age spread (Baraffe et al. 2009, 2012), and can account for unexpected lithium and beryllium depletion in some young stars (Baraffe & Chabrier 2010; Viallet & Baraffe 2012). In addition, quiescent periods between the accretion bursts can promote disk fragmentation, which might in turn promote planet formation (Stamatellos et al. 2011).

Episodic accretion may also have important implications for the chemical evolution of protostellar disks and envelopes. Using simplified core collapse calculations with prescribed accretion bursts, Lee (2007) and Visser & Bergin (2012) showed that the bursts can lead to evaporation of CO, CO2, and some other ices in protostellar envelopes. The freeze-out timescale of these species is expected to be longer than the typical duration of the burst (~100 yr), so that the chemical signatures of the bursts can linger through the quiescent phase of accretion. Kim et al. (2011, 2012) found evidence for pure CO2 ice in about 50% of the observed low-luminosity sources, implying that the dust temperature must have been higher in the past presumably due to recent accretion bursts.

In this paper, we employ numerical hydrodynamics simulations coupled with a stellar evolution code, as done in Baraffe et al. (2012), and a simplified chemical model, to explore the effect of episodic accretion onto the evaporation and freeze-out of CO and CO2 during the early stages of star formation. The coupling between hydrodynamical simulations and the evolution of the central, accreting object provides a consistent value of the protostar luminosity, which is relevant for the estimate of the radiative feedback from the protostar in the disk and the envelope. This feedback impacts the gas/dust temperatures which are required for the calculations of CO and CO2 abundances. In the present hydrodynamical calculations, estimate of the gas/dust temperatures relies on the diffusion approximation for the treatment of radiative transfer and additional approximations regarding the geometry of the disk and of the surrounding material. We analyse the uncertainties of temperatures based on such approximations, by comparing them with temperatures obtained from improved radiative transfer calculations based on the TORUS code (Harries et al. 2004; Harries 2011). We will also discuss (Sect. 3) the impact of our approximation on the predicted CO and CO2 abundances.

2. Model description

2.1. Basic equations

Our numerical hydrodynamics model for the formation and evolution of a young stellar object is described in detail in Vorobyov & Basu (2010). Here, we briefly review the main concepts and describe latest modifications. The model in the most general case includes a forming protostar, described by a stellar evolution code, and a protostellar disk plus infalling envelope, both described by numerical hydrodynamics equations. We use the thin-disk approximation, which is an excellent means to calculate the evolution for many orbital periods and many model parameters and its justification is discussed in Vorobyov & Basu (2010). The thin-disk approximation is complemented by a calculation of the vertical scale height h in both the disk and envelope determined in each computational cell using an assumption of local hydrostatic equilibrium. The resulting model has a flared structure with the vertical scale height increasing with radial distance according to the law h ∝ r1.5. Both the disk and envelope receive a fraction of the irradiation energy from the central protostar described by Eq. (7) (see below). The main physical processes taken into account when computing the evolution of the disk and envelope include viscous and shock heating, irradiation by the forming star, background irradiation, radiative cooling from the disk surface and self-gravity. The corresponding equations of mass, momentum, and energy transport are

where subscripts p and p′ refers to the planar components (r,φ) in polar coordinates, Σ is the mass surface density, e is the internal energy per surface area, is the vertically integrated gas pressure calculated via the ideal equation of state as with γ = 7/5, is the velocity in the disk plane, and is the gradient along the planar coordinates of the disk. The gravitational acceleration in the disk plane, , takes self-gravity of the disk into account, found by solving for the Poisson integral (see details in Vorobyov & Basu 2010), and the gravity of the central protostar when formed.

Turbulent viscosity due to sources other than gravity is taken into account via the viscous stress tensor Π, the expression for which is provided in Vorobyov & Basu (2010). We parameterize the magnitude of kinematic viscosity ν using the α-prescription with a spatially and temporally uniform α = 5 × 10-3.

The radiative cooling Λ in Eq. (2) is determined using the diffusion approximation of the vertical radiation transport in a one-zone model of the vertical disk structure (Johnson & Gammie 2003) (4)where σ is the Stefan-Boltzmann constant, is the midplane temperature of gas1, μ = 2.33 is the mean molecular weight, R is the universal gas constant, and ℱc = 2 + 20tan-1(τ)/(3π) is a function that secures a correct transition between the optically thick and optically thin regimes. We use frequency-integrated opacities of Bell & Lin (1991). The heating function is expressed as (5)where Tirr is the irradiation temperature at the disk surface determined by the stellar and background black-body irradiation as (6)where Tbg is the uniform background temperature (in our model set to the initial temperature of the natal cloud core) and Firr(r) is the radiation flux (energy per unit time per unit surface area) absorbed by the disk surface at radial distance r from the central star. The latter quantity is calculated as (7)where γirr is the incidence angle of radiation arriving at the disk surface at radial distance r. The incidence angle is calculated using the disk surface curvature inferred from the radial profile of the disk vertical scale height (see Vorobyov & Basu 2010, for more details).

The central object’s luminosity L∗ is the sum of the accretion luminosity L∗ ,accr = GM∗Ṁ/2R∗ arising from the gravitational energy of accreted gas and the photospheric luminosity L∗ ,ph due to gravitational contraction and deuterium burning in the protostar interior. The stellar mass M∗ and accretion rate onto the star Ṁ are determined self-consistently during numerical simulations using the amount of gas passing through the sink cell. The evolution of the accreting protostar is based on the Lyon stellar evolution code with input physics described in Chabrier & Baraffe (1997) and including accretion processes as described in Baraffe et al. (2009, 2012). The accretion rates are derived from the hydrodynamic calculations above described. As in Baraffe et al. (2012), we assume that a fraction α of the accretion energy is absorbed by the protostar, while the fraction (1-α) is radiated away and contributes to L∗ in Eq. (7)2. In the present calculations, we adopt a “hybrid” scheme to describe the contribution of the accreted matter to the protostar’s internal energy, with “cold” accretion, i.e. α = 0, when accretion rates remain smaller than a critical value Ṁcr, and “hot” accretion, i.e. α ≠ 0, when Ṁ > Ṁcr. In this paper we adopt Ṁcr = 10-5  M⊙ yr-1 and α = 0.2 (see discussion in Baraffe et al. 2012). For the initial mass of the protostar, corresponding to the second Larson core mass, we adopt a value of 1.0 MJup with an initial radius ~1.0  R⊙.

The input parameters provided by the hydrodynamic calculations to the coupled stellar evolution code are the age of the protostar and the accretion rate onto the stellar surface Ṁ, while the output are the radius and the photospheric luminosity of the protostar. The stellar evolution code is called to update the properties of the protostar every 5 yr of the physical time. For comparison, the global hydrodynamical timestep may be as small as a few weeks and the whole cycle of numerical simulations may exceed 1.0 Myr. This coupling of the disk and protostar evolution allows for a self-consistent determination of the radiative input of the protostar into the disk thermal balance, which is important for the accurate study of disk instability and fragmentation.

The dynamics of carbon monoxide (CO) and carbon dioxide (CO2) including absorption onto and desorption from dust grains is computed using the modified equations of continuity for the surface densities of solid () and gaseous () phases where index i corresponds to either CO or CO2, λ is the absorption rate (s-1) and η is the desorption velocity from dust grains in units of g cm-2 s-1. The expressions for η and λ were taken from Charnley et al. (2001) and Visser et al. (2009) assuming a zero-order desorption for thick ice mantles

where Tmp is the gas midplane temperature3, Mi and mi are the molecular weight and mass (in gram), respectively, and nd is the number density of dust grains. The values of the binding energy for CO and CO2 (Eco/k = 885 and Eco2/k = 2300 K) and vibrational frequency (νco = 7 × 1026 cm-2 s-1 and νco2 = 9.5 × 1026 cm-2 s-1) were taken from Bisschop et al. (2006) and Noble et al. (2012). In this paper, we make a simplifying assumption of Td = Tmp and discuss possible consequences in Sect. 4.

In the limit of a constant radius of dust grains a, the expression in brackets in Eq. (10) can be written in the following simplified form (12)where Ad2g = 0.01 is the dust to gas mass ratio, md = 4πa3ρd.p./3 is the mass of a dust grain and ρd.p. = 2.5 g cm-3 is the density of dust grains. The radius of dust grains is set in this study to a = 0.1  μm.

When writing Eqs. (8) and (9) for the surface densities of CO and CO2 we assumed that dust is passively transported with gas and neglected the dust radial drift caused by the dust–gas drag force. The latter simplification is of little consequence for the dynamics of dust grains with sizes ≲10 μm on timescales of ≤0.5 Myr (Takeuchi & Lin 2002). We also do not take chemical reactions that may change the abundance of CO and CO2 into account, and consider only phase transformations of these species. Chemical transformations will be taken into account in a follow-up study.

We start our numerical simulations from the gravitational collapse of a starless cloud core, continue into the embedded phase of star formation, during which a star, disk, and envelope are formed, and terminate our simulations in the T Tauri phase, when most of the envelope has accreted onto the forming star plus disk system. The protostellar disk, when present, is located in the inner part of the numerical polar grid, while the collapsing envelope occupies the rest of the grid. As a result, the disk is not isolated but is exposed to intense mass loading from the envelope. In addition, the mass accretion rate onto the disk is not a free parameter of the model but is self-consistently determined by the gas dynamics in the envelope.

To avoid too small time steps, we introduce a “sink cell” at rsc = 5 AU and impose a free inflow inner boundary condition and free outflow outer boundary condition so that the matter is allowed to flow out of the computational domain but is prevented from flowing in. The sink cell is dynamically inactive; it contributes only to the total gravitational potential and secures a smooth behaviour of the gravity force down to the stellar surface. During the early stages of the core collapse, we monitor the gas surface density in the sink cell and when its value exceeds a critical value for the transition from isothermal to adiabatic evolution, we introduce a central point-mass object. In the subsequent evolution, 90% of the gas that crosses the inner boundary is assumed to land onto the protostar. A small fraction of this mass (a few per cent) remains in the sink cell to guarantee a smooth transition of the gas surface density across the inner boundary. The other 10% of the accreted gas is assumed to be carried away with protostellar jets.

The numerical resolution is 512 × 512 grid points and the numerical procedure to solve hydrodynamics Eqs. (1)−(2) is described in detail in Vorobyov & Basu (2010). The modified continuity Eqs. (8) and (9) are solved using a two-step procedure. First, the advection part (zero r.h.s.) is solved using the same piecewise parabolic advection scheme as for the gas surface density. Then, the surface densities of CO and CO2 are updated to take the phase transformations into account using a first-order backward Euler scheme. This scheme is implicit and in principle does not require a time-step limiter. However, too big timesteps can result in the loss of accuracy. Therefore, if the relative change of and/or exceeds 10% over one global hydro step, the local integration timestep for Eqs. (8) and (9) is reduced by a factor of 2 and the solution is sought once again. This subcycling procedure is repeated until the desired accuracy is achieved.

2.2. Initial setup

For the gas surface density Σ and angular velocity Ω radial distributions we take those typical of pre-stellar cores formed as a result of the slow expulsion of magnetic field due to ambipolar diffusion, with the angular momentum remaining constant during axially-symmetric core compression (Basu 1997)

Here, Ω0 and Σ0 are the angular velocity and gas surface density at the centre of the core and is the radius of the central plateau, where cs is the initial sound speed in the core. The gas surface density distribution described by Eq. (14) can be obtained (to within a factor of order unity) by integrating the three-dimensional gas density distribution characteristic of Bonnor-Ebert spheres with a positive density-perturbation amplitude A (Dapp & Basu 2009). The value of A is set to 1.2 and the initial gas temperature is set to 10 K.

In order to form a gravitationally unstable core, we set the ratio of the outer radius rout to the radius of the central plateau r0 to 6.0. For the outer radius of the core rout = 16   000 AU, the resulting central surface density is Σ0 = 4.5 × 10-2 g cm-2 and the total mass of the core is Mc = 1.23  M⊙. The central angular velocity Ω0 is set to 1.0 km s-1 pc-1, which yields the ratio of rotational to gravitational energy β = 5 × 10-3. Our previous numerical simulations indicate that pre-stellar cores with similar Mc and β produce disks that are gravitationally unstable. Mass accretion rates onto the protostar in these models often exhibit episodic accretion bursts caused by disk fragmentation and migration of the fragments onto the protostar (Vorobyov & Basu 2010).

The abundances of CO and CO2 (relative to the number density of molecular hydrogen) are set to 5 × 10-5, typical for pre-stellar gravitationally unstable cores4 (Charnley et al. 2001; Pontoppidan et al. 2008). We assume that initially 10% of CO and CO2 is in the solid phase in order to be consistent with previous works of Visser et al. (2009) and Visser & Bergin (2012).

thumbnail Fig. 1

Gas surface density images in the inner 1000 × 1000 AU showing the disk evolution during 0.24 Myr after the formation of the central star. Vigorous disk fragmentation is evident in the figure. The scale bar is in log g cm-2.

3. Results of numerical simulations

3.1. Main characteristics of the model

We start by describing the main properties of the forming protostar, disk and envelope. Figure 1 shows the gas surface density maps in the inner 1000 × 1000 AU box at six consecutive times after formation of the central protostar. The disk forms at t = 0.01 Myr, grows in mass and size due to continuing mass loading from the infalling envelope, and becomes gravitationally unstable (as manifested by a weak spiral structure) as early as at t = 0.04 Myr. The first episode of disk fragmentation occurs at t ≈ 0.08 Myr. In the subsequent evolution, multiple fragments emerge in the disk at distances ≳50 AU but most migrate into the inner regions and through the sink cell due to the loss of angular momentum via gravitational interaction with the spiral arms and other fragments in the disk. This phenomenon is studied in detail by Vorobyov & Basu (2006, 2010) and is confirmed by independent fully three-dimensional studies (e.g. Machida et al. 2011a). The ultimate fate of the migrating fragments depends on how quickly they can contract to stellar or planetary-sized objects. If the contraction timescale is longer than the migration/tidal destruction timescale, then the fragments will be completely destroyed when approaching the protostar, releasing their gravitational energy in the form of luminosity outbursts comparable to magnitude to FU Orionis and EX Lupi objects. This scenario is assumed in the present work.

Figure 2 presents the time evolution of a) the mass accretion rate onto the star Ṁ; b) the photospheric (dashed line) and total (solid line) luminosities; c) the masses of the protostar (solid line), protostellar disk (dashed line), and envelope (dash-dotted line); and d) the radius of the disk. The partition between the disk and infalling envelope is based on a threshold density of Σd2e = 0.5 g cm-2 and the corresponding algorithm is described in detail in Dunham & Vorobyov (2012). The time is counted from the formation of the central protostar.

Two luminosity outbursts exceeding in magnitude 100  L⊙, along with a number of smaller bursts (≳several × 10  L⊙), are evident against the background luminosity of a few L⊙. The photospheric luminosity provides a negligible input into the total flux in the early evolution (t < 0.04 Myr) because the accretion rate does not exceed 10-5  M⊙ yr-1 and the accretion process is essentially cold, depositing little entropy to the protostar The first episode of hot accretion occurs at t ≈ 0.04 Myr, when the mass accretion rate exceeds 10-5  M⊙ yr -1. Thereafter, the photospheric luminosity is slowly varying around 2 − 3  L⊙.

thumbnail Fig. 2

Main characteristics of the forming system as a function of time elapsed since the formation of the central star: a) the mass accretion rate onto the protostar; b) the total and photospehric luminosities; c) the masses of the protostar, disk, and envelope; and d) the disk radius.

Panels c) and d) in Fig. 2 illustrate the time evolution of other global characteristics of our model. The disk mass steadily grows during the initial evolution and reaches a value of Md ≈ 0.35  M⊙ at 0.2 Myr. In the subsequent evolution, the disk mass saturates and starts to decline slowly. The mass of the protostar exceeds that of the disk and shows episodic sharp increases caused by massive fragments spiralling in onto the protostar. The envelope mass steadily decreases and the model enters the class I stage of stellar evolution, defined as the time when approximately half of the initial mass reservoir is left in the envelope, at t ≈ 0.11 Myr. The disk radius increases with time but also shows significant radial pulsations. These contractions/expansions are caused by migration of the fragments and the corresponding redistribution of angular momentum within the disk, when the fragment’s angular momentum is transferred to the spiral arms causing them to unwind and expand radially outward. Overall, the considered model produces a rather massive and extended disk. Theoretically, such disks are not unexpected in the embedded stage of star formation characterized by high rates of mass infall onto the young disk. However, it must be kept in mind that we have not taken into account magnetic fields and the effect of the external environment, which can significantly decrease the mass and size of protostellar disks (see e.g. Hennebelle & Teyssier 2008). Therefore, although the qualitative picture issued from the present calculations should be rather robust, the quantitative results must be considered with due caution.

thumbnail Fig. 3

Effect of the luminosity bursts on the spatial distribution of gas-phase CO and CO2. The four columns present (from left to right) the gas surface density distribution (log g cm-2), gas temperature (K), CO gas-phase fraction, and CO2 gas-phase fraction in the inner 3000 × 3000 AU region. The four rows show (from top to bottom) the pre-burst phase (t = 0.128 Myr), the burst phase (t = 0.13 Myr), the post-burst phase (t = 0.134 Myr) and the quiescent phase (t = 0.14 Myr). The black lines are isotemperature contours outlining the regions where gas temperature exceeds 20 K (third column) and 40 K (fourth column) above which CO and CO2 are supposed to be in the gas phase in an equilibrium, time-independent case.

3.2. Spatial distribution of CO and CO2

Figure 3 presents (from left to right) the gas surface density Σ (first column), the gas temperature Tmp (second column), the CO gas-phase fraction (third column) and the CO2 gas-phase fraction (fourth column) in the inner 3000 × 3000 AU region, at four distinct time instances. The first row (from top to bottom) highlights the evolution stage soon after a moderate luminosity burst with L∗ ≈ 20  L⊙ (Ṁ ≈ 10-5  M⊙ yr-1) that occurred at t = 0.127 Myr. The second row presents the model during a strong luminosity burst with L∗ ≈ 250  L⊙ (Ṁ ≈ 2 × 10-4  M⊙ yr-1) at t = 0.13 Myr. The third row shows the model at t = 0.134 Myr, i.e., 4.0 kyr after the strong burst. Finally, the fourth row represents the quiescent stage at t = 0.14 Myr with L∗ ≈ 3.5  L⊙. The black lines delineate isotemperature contours with Tmp = 20 K (third column) and Tmp = 40 K (fourth column). These values correspond to the evaporation temperatures of CO and CO2 ices from dust grains (Noble et al. 2012).

thumbnail Fig. 4

CO gas-phase fraction (red lines) and total stellar luminosity L∗ (black lines) vs. time. Two time periods (20 kyr each) are shown to highlight two strongest luminosity bursts and several moderate ones. The correlation between and L∗ is evident. In particular, steeply rises during the burst to a maximum value and gradually declines to a minimum value after the burst.

A visual inspection of Fig. 3 indicates that the protostellar disk (localized in the inner 300−350 AU) is characterized by almost complete evaporation of CO from dust grains (), a result consistent with the recent study of molecular abundances in gravitationally unstable disks by Ilee et al. (2011). This holds for both the burst and quiescent stages of disk evolution. A moderate burst at t = 0.127 Myr (L∗ ≈ 20  L∗) can also evaporate CO in the innermost parts of the envelope up to r ≈ 500 AU (first row) and the CO evaporation region can extend beyond 1000 AU for a strong burst with L∗ ≈ 250  L⊙ (second row).

The burst duration is usually limited to 100−200 yr (defined as the period of time during which the mass accretion rate constantly exceeds 10-5  M⊙ yr-1) but the effect of the burst is lingering in the envelope for a significantly longer time. For instance, the duration of the burst at t = 0.13 Myr is about 0.2 kyr, but a significant amount of gas-phase CO is seen in the envelope at t = 0.134 Myr, i.e., almost 4.0 kyr after the burst (third row in Fig. 3). The gas temperature in the envelope at this stage drops below 20 K (as indicated by the isotemperature contour), which means that the presence of gas-phase CO is truly related to a recent luminosity burst. Two lobe-like features that are almost completely devoid of gas-phase CO are evident in the third row. The spiral arms have dynamically swept through the lobe regions prior to the snapshot, and during this passage the CO gas has frozen out onto the grains in the high-density wake of the arms. At the same time, the leading edge of the arms is shock compressed to a sufficiently high temperature to evaporate CO. Finally, we note that CO freezes out onto dust grains during the quiescent stage if its duration is comparable to or longer than 10 kyr (fourth row in Fig. 3).

The effect of a recent luminosity burst on the gas-phase CO in the envelope can be understood from the following simple analysis. The e-folding time of freeze-out onto dust grains tads = λ-1 can be expressed as (15)where ρg = Σ/2h is the gas volume density. For the typical conditions in the inner envelope at r = 1000 AU (Tmp = 15 K and ng = 106 cm-3), the resulting e-folding time for CO freeze-out onto dust grains is approximately 2500 yr. The fact that tads for carbon monoxide can be much longer than the burst duration opens up a possibility for the observational detection of recent bursts, as suggested by Lee (2007) and Visser & Bergin (2012) on the basis of simplified core collapse calculations. Equation (15) also indicates that tads increases linearly with the radius of dust grains a, suggesting that this phenomenon may become even more pronounced if dust grains have enough time to coagulate and grow to sizes greater than adopted in the present study, a = 0.1  μm.

In contrast, the phase transformations of CO2 during the bursts are much less pronounced. Most of the disk and all of the envelope in Fig. 3 are characterized by CO2 frozen out onto grains. CO2 has an evaporation temperature of 35−40 K and in the quiescent stage the gas-phase CO2 is present only in the inner 25−30 AU. During the strong burst with L∗ ≈ 250  L⊙ (second row in Fig. 3), the temperature may rise above the CO2 evaporation limit in the inner 200−300 AU, transforming most of CO2 into the gas phase. However, the gas density in this region of the disk (ng = (1 − 5) × 109 cm-3) is considerably higher than in the envelope and the e-folding time for the CO2 freeze-out onto dust grains (tads ~ 1.0 yr) is significantly shorter than that of carbon monoxide (tads ~ 2500 yr). As a result, the gas-phase CO2 quickly returns into the solid phase after the burst and no gas-phase CO2 is seen at t = 0.134 Myr at radial distances beyond 25−30 AU. We conclude that the phase transitions of CO2, and in particular the abundance of gas-phase CO2, are less convenient for monitoring the recent burst activity. We note, however, that the abundance of solid CO2 appears to be sensitive to the past accretion history and can be a good episodic accretion tracer, as recently demonstrated by Kim et al. (2011).

3.3. Time variations of the gas-phase CO and CO2 in the envelope

In Fig. 4 we consider the time evolution of the gas-phase CO fraction () during two time intervals: 0.08−0.1 Myr (left column) and 0.12−0.14 Myr (right column). These time intervals were chosen so as to capture two most intense and several moderate luminosity bursts. The red lines present averaged over three radial bins: 300 < r < 1000 AU (top row), 300 < r < 2000 AU (middle row), and 500 < r < 2000 AU (bottom row). The first two bins cover the outer parts of the disk and the inner parts of the infalling envelope, while the last bin covers only the inner envelope (the disk radius is less than 500 AU, see Fig. 2). The time is counted from the formation of the central protostar. The black lines show the total protostellar luminosity L∗ vs. time.

The general correlation between L∗ and is evident in Fig. 4. The CO gas-phase fraction steeply rises during the burst to a maximum value and gradually declines to a minimum value after the burst. The relaxation time to the pre-burst stage is notably longer than the burst duration, in agreement with analytic estimates performed in the previous section. This pattern of behaviour – periodic steep rises during the bursts followed by gradual declines – is most pronounced in the envelope. Carbon monoxide in the disk is mostly in the gas phase and the bursts have little effect on the CO gas-phase fraction there. The top-right panel in Fig. 4 illustrates such an example. The CO gas-phase fraction is averaged over the 300 < r < 1000 AU radial bin, which covers part of the disk (with radius 350−400 AU at this time instant). Evidently, the correlation between and L∗ is less pronounced. It is therefore important to filter out the disk contribution when calculating .

thumbnail Fig. 5

Top: CO gas-phase fraction vs. gas temperature Tmp calculated in the 500 < r < 1000 AU radial annulus during a time period covering a strong luminosity burst at t = 0.13 Myr. The dashed, solid, and dash-dotted lines show the pre-burst, burst, and post-burst phases with durations of 1.0 kyr, 0.2 kyr, and 5 kyr, respectively. The arrows indicate this evolution sequence. The horizontal and vertical dotted lines are the CO fraction of 0.5 and the gas temperature at which the CO desorbs. The grey-shaded area highlights the phase-space region with an abnormally high fraction of the gas-phase CO. The thick dash-dotted line highlights the –Tmp model track in this region. Bottom: CO2 gas-phase fraction vs. Tmp calculated in the 100 < r < 300 AU radial annulus during a strong luminosity burst at t = 0.082 Myr. The dashed, solid, and dash-dotted lines show the pre-burst, burst, and post-burst phases with durations of 0.4 kyr, 0.2 kyr, and 0.2 kyr, respectively. The horizontal and vertical dotted lines mark the CO2 fraction of 0.5 and the gas temperature at which the CO2 desorbs. The grey-shaded area highlights the phase-space region with an abnormally high fraction of the gas-phase CO2.

The delayed adsorption of CO onto dust grains after a strong luminosity burst is illustrated in the top panel of Fig. 5 showing the gas-phase CO fraction vs. gas temperature Tmp in the 500 < r < 1000 AU radial annulus. A short time period covering the burst at t = 0.13 Myr is shown including 1.0 kyr immediately before the burst (dashed line), 0.2 kyr during the burst (solid line), and 5 kyr after the burst (dash-dotted line). The vertical and horizontal dotted lines mark the CO fraction of 0.5 and the gas temperature above which CO evaporates. It is evident that the evolution tracks of the CO gas-phase fraction are different in the pre- and post-burst phases. Before the burst CO is found mostly in the solid phase, as expected for Tmp < 20 K. During the burst, increases with the rising gas temperature and reaches a maximum value of ≈1.0. However, after the burst the gas temperature quickly drops below 20 K, but remains abnormally high. The corresponding region in the –Tmp diagram is filled with a grey shade and the model track in this region is highlighted by the thick dash-dotted line. The time spent in this abnormal region (with and Tmp < 20 K) is 1.2 kyr5 The resulting mismatch between the gas temperature and the CO gas-phase content in the post-burst phase can be used to detect recent luminosity outbursts. Objects that are found in the grey-shaded area of the diagram are likely to be in the post-burst phase; objects in the quiescent and pre-burst phases are unlikely to fall into this region.

On the contrary, the phase transformation of CO2 lacks such a characteristic feature. The bottom panel in Fig. 5 presents the gas-phase CO2 fraction vs. gas temperature Tmp in the 100 < r < 300 AU radial annulus during a short time period covering the burst at t = 0.082 Myr. The dotted vertical line marks the critical gas temperature above/below which CO2 is supposed to be mostly in the gas/solid phase. During the burst rises to a maximum value of ≈1.0 and returns to a small value of ≤0.1 during just 0.2 kyr in the post-burst phase. There is no evolution stage when the gas temperature is below the condensation temperature of CO2, but the gas-phase CO2 fraction is abnormally high. The evolution tracks of the CO2 gas-phase fraction in the pre- and post-burst phases follow a similar path in the Tmp– diagram and do not pass through the shaded region. This makes CO2 ineffective in detecting recent luminosity bursts.

4. Model uncertainties

Our model has a few assumptions that may affect the abundance of CO and CO2 ices in the disk and envelope. In this section, we discuss the model caveats and their impact on our main results.

Uncertainties in the gas/dust temperature. In our model we assumed that the gas and dust temperatures are equal. This is usually correct for the gas densities typical for a protostellar disk (>109 − 1010 cm-3), where frequent collisions between the gas and dust particles lead to fast thermalization between these species. In the envelope, however, the gas temperature may differ from that of the dust. A more accurate approach with a separate treatment of dust and gas thermodynamics is needed to assess the possible effect of the gas to dust temperature imbalance.

The lack of vertical structure. The use of the thin-disk approximation implies that all quantities are averaged over the scale height and there are no variations in the gas volume density and temperature with the distance from the midplane. This may affect the adsorption and desorption rates of ices if strong vertical variations in the thermodynamical properties are present. Work is in progress to reconstruct the vertical structure based on the coupled solution of the radiation transfer and hydrostatic equilibrium equations. Preliminary results indicate that notable vertical variations in the gas temperature and density may be present in the inner disk regions but they are diminishing in the envelope (Vorobyov et al., in prep.).

Simplified treatment of gas thermodynamics. Present results are based on a simplified form of diffusion approximation to calculate the thermal balance in the disk and envelope. We test our approach by comparing our derived gas/dust temperatures with those obtained with an improved treatment of radiative transfer using the TORUS code (Harries et al. 2004; Harries 2011). The later employs the Monte-Carlo algorithm described by Lucy (1999). We assume silicate dust grains with a dust-to-gas ratio of 0.01, and adopt a canonical ISM grain-size distribution (Mathis et al. 1977). The grain opacities and Mie phase matrices were calculated from the refractive indices of astronomical silicates (Draine & Lee 1984). Note that TORUS also assumes that dust and gas temperatures are the same.

thumbnail Fig. 6

Top: gas radial temperature profiles. In particular, the red line presents the azimuthally averaged midplane temperature Tmp derived from hydrodynamical simulations at t = 0.131 Myr and the total stellar luminosity L∗ = 9  L⊙. The black solid line shows the gas temperature derived using the TORUS radiation transfer code, for the same gas density distribution and stellar luminosity as in the hydro model (case 1). The dashed solid line presents the TORUS calculation of the gas temperature for the same gas density distribution (case 1) but a higher stellar luminosity of 100  L⊙. Finally, the blue line shows the gas temperature calculated by the TORUS code for a spherically symmetric gas density distribution representing an infalling spherical envelope with the same mass as in the hydro model (case 2). Bottom: the gas volume density distribution used by the TORUS code (case 1) to calculate the gas temperatures (black solid and dashed lines). The scale bar is in g cm-3.

The post-processed radiative transfer calculations with TORUS require the knowledge of the density distribution. We explore two different cases with axisymmetric (case 1) and spherically symmetric (case 2) distributions, respectively. Case 1 provides information on the temperature uncertainty resulting from the simplified treatment of radiative transfer in the hydrodynamical model. Case 2 corresponds to a spherically symmetric envelope and provides a limiting case which helps estimating the temperature uncertainty in the envelope due to imposed thin-disk limit (see below and Fig. 6). The combination of both cases can thus provides a rough estimate of the error on temperatures determined in the hydrodynamical model.

The top panel in Fig. 6 presents the gas temperature vs. radial distance calculated using both our hydrodynamical model (providing the so-called midplane gas temperature Tmp) and the TORUS code. In particular, the red line shows the azimuthally averaged Tmp in our model in the post-burst stage at t = 0.131 Myr when the luminosity of the central protostar has dropped to approximately 9.0 L⊙. The solid black line (case 1) is the gas midplane temperature calculated using the TORUS code for the same gas density distribution and the same stellar luminosity as in our hydrodynamical model at t = 0.131 Myr. A mass-weighted vertical mean temperature differs from the midplane temperature in the TORUS calculation by a factor of order unity. We convert the model surface densities into volume densities (needed in TORUS) using the following simple formula: , where and are the azimuthally averaged gas surface density and vertical scale height. The resulting gas density distribution ρ(r) has a flared form as shown in the bottom panel of Fig. 6, with the vertical scale height increasing with distance. The color shaded area represents the disk plus envelope, while the white area is filled with a rarefied gas representing an outflow cavity.

The dashed black line (case 1) in top panel of Fig. 6 corresponds to the TORUS calculation of the gas temperature for the same model gas density distribution but for a higher luminosity of the central source during the burst, namely L∗ = 100  L⊙. A comparison between those different temperatures is mostly relevant for the envelope (i.e r ≳ 300 AU), since the radiative equilibrium calculations with TORUS only include heating from the central protostar and no dynamical effects of the gas (e.g. viscous heating or compression), which are more important in the disk than in the envelope. This explains the significantly lower temperature obtained in the inner region (r ≲ 300 AU) with TORUS compared to Tmp, for similar gas density distribution and protostar luminosity. The coupling of the hydrodynamical model with TORUS, providing a more consistent description of heating and cooling processes, is left for future work.

In the outer regions (r ≳ 300 AU) the temperature differences between the TORUS temperature (case 1, solid black line) and Tmp is less than 30%. Though small, such temperature uncertainty may be relevant if temperatures are close to condensation temperatures, as indicated by the horizontal dotted line for CO condensation in the upper panel of Fig. 6. Reassuringly enough, case 1 with a higher protostar luminosity (dashed black line) shows the overall temperature increase expected during a burst. Such a temperature increase will push regions where CO may condense much further out (≫1000 AU), giving some confidence in the abundance estimate of the gas-phase CO in the envelope derived in Sect. 2.3. Moreover, an overall decrease in the gas density with radial distance implies hat the typical freeze-out time of gas-phase CO onto dust grains will increase at larger distances, enabling the detection of earlier bursts.

Finally, the blue line in the top panel of Fig. 6 shows the gas temperature calculated by the TORUS code for a spherically symmetric gas density distribution (case 2) representing an infalling spherical envelope with mass equal to the envelope mass in our hydrodynamic model (0.5  M⊙). The gas density profile in the envelope follows the relation ρ(r) ∝ r-2, characteristic for collapsing truncated cores of finite size (Vorobyov & Basu 2005). The inner and outer radii of the envelope are set to 100 AU and 15 000 AU. The temperature difference resulting from different envelope geometries (TORUS case 1 versus TORUS case 2) can be large at ~100 AU (about a factor three), but rapidly decreases as a function of radial distance. A spherical envelope is an extreme case, since outflows are expected to create a cavity which can extend as far as several thousand AU or more (Commerçon et al. 2010; Machida et al. 2011b). The spherical case is however interesting since it provides an upper limit to the expected temperature in the envelope. Current assumptions made in the hydrodynamical model are thus not expected to alter the conclusions derived from the predicted gas-phase abundances of CO during and after a burst (Sect. 2.3). But the spherical case shows that taking the vertical structure of the disk and the envelope into account is important for a better estimate of the temperature.

The lack of chemical reactions. In the present work, we neglect reactions that can alter the chemical composition of protostellar disks and collapsing envelopes. In particular, Kim et al. (2011) suggest that CO turns into CO2 during the quiescent accretion stage when CO freezes out onto dust grains. They found that up to 80% of the original CO content can be converted into CO2 in the envelope during the main accretion phase. We note, however, that the abundance of the gas-phase CO in the quiescent stage is non-negligible (see Fig. 4), owing to a long characteristic freeze-out time, and this may reduce the efficiency of CO to CO2 conversion on the surface of dust grains. Detailed calculations using the grain-surface reactions are needed to estimate the CO depletion due to this effect.

The hybrid cold/hot accretion scheme. In this study, we assumed that a fraction of the accretion energy is absorbed by the protostar when the accretion rate Ṁ exceeds a critical value of Ṁcr = 10-5  M⊙ yr-1. Below this value, accretion onto the protostar proceeds in the cold regime, with all accretion energy radiated away. The choice of Ṁcr was motivated by fitting of the spectral energy distribution of FU Orionis by Hartmann et al. (2011), who argued that these eruptive stars undergo expansion during the burst, which in turn indicates that some accretion energy is absorbed by the star. Accretion rates during FU-Orionis-type bursts typically lie in the 10-6 − 10-4  M⊙ yr-1 range. Varying the value of Ṁcr in these limits, will affect the stellar radius R∗ and the resulting accretion and photospheric luminosities. However, the total stellar luminosity is not expected to vary significantly because the increase in the accretion luminosity due to, say, decrease in the stellar radius, will be at least partly compensated by the corresponding decrease in the photospheric luminosity.

5. Conclusions

In this paper, we studied the phase transformations of CO and CO2 during the early stages of low-mass star formation characterized by strong luminosity outbursts resulting from accretion bursts. We used numerical hydrodynamics modelling in the thin-disk limit to describe the gravitational collapse of a rotating pre-stellar core, extending our calculations into the disk formation stage and terminating the calculations when most of the initial pre-stellar core has accreted onto the disk plus protostar system. In the early accretion phase, the system experiences repetitive accretion and luminosity bursts caused by disk gravitational fragmentation and quick inward migration of the fragments onto the protostar (Vorobyov & Basu 2010). The basic hydrodynamics equations were complemented with the continuity equations describing the adsorption and desorption of CO and CO2 onto and from dust grains. We calculated the gas-phase fractions of CO and CO2 ( and ) in the protostellar disk and infalling envelope in the pre-burst, burst, and post-burst phases. We found the following.

  • In the quiescent phase characterized by total stellar luminosity ofthe order of a few L⊙, CO in the disk is found mostly in the gas phase, while in the envelope CO has mostly frozen out onto dust grains. CO2 is found in the gas-phase only in the inner 20–30 AU, while in the rest of the disk and envelope CO2 freezes out onto dust.

  • During strong luminosity bursts characterized by a total luminosity from a few tens to a few hundreds L⊙, CO ice evaporates from dust grains in part of the envelope, while CO2 ice turns into the gas phase only in the inner several hundred AU of the disk.

  • The typical time for freeze-out of the gas-phase CO onto dust grains in the envelope (a few kyr) is considerably longer than the typical duration of a luminosity burst (0.1–0.2 kyr). As a result, a significant amount of the gas-phase CO can still be present in the envelope long after the system has returned into the quiescent phase. This phenomenon can be used to infer recent luminosity bursts with magnitudes typical of EX-Lupi and FU-Orionis-type outbursts, as suggested by recent semi-analytical studies by Lee (2007) and Visser & Bergin (2012).

  • In contrast, the typical freeze-out time of the gas-phase CO2 is comparable to the burst duration, owing to significantly higher gas densities in the disk. We thus conclude that CO2 is probably not a good tracer for recent burst activity in young protostars.

Regarding uncertainties of the present model, the heating due to stellar irradiation is crucial to derive temperatures in the disk and the envelope. In this work, we tried to estimate uncertainties resulting from approximate treatment of radiative transfer in the hydrodynamical model by comparing temperatures obtained from an improved treatment of irradiation effects based on a radiative transfer code. The results of this comparison indicate that the uncertainty on the predicted temperatures depends on the radial distance from the protostar and the structure of the envelope. However, this uncertainty does not alter our conclusions regarding the abundance of gas-phase CO in the envelope and the possibility to use it as a tracer of recent accretion burst activity. Our short-term plans to reduce the present uncertainties involving (i) the improvement of radiative transfer calculations, based on the coupling between the hydrodynamical model and a radiative transfer code, like e.g. TORUS; (ii) the reconstruction of the vertical structure of the disk, in order to improve over the thin disk approximation and (iii) the implementation of a chemical model, which is more challenging but would provide more robust predictions. Finally, an important weakness of our model is the absence of magnetic fields, which are expected to alter the properties of the collapse, the disk (e.g. the size) and the accretion process onto the protostar. Inclusion of magnetic field, with the improvement above-mentioned, is a major challenge that will provide the state-of-the-art on a longer term. Three-dimensional non-ideal MHD simulations including radiation hydrodynamics are now underway to explore the second collapse (Masson et al. 2012). The work presented here already provides consistent predictions that can be tested against observations. The increased sophistication on the MHD, thermal and chemical treatments we propose will further improve our understanding of embedded phases and episodic accretion.

1

This definition of the midplane temperature is accurate within a factor of order unity (Zhu et al. 2012)

2

As in Baraffe et al. (2009), we assume a value ϵ = 1/2 characteristic of accretion from a thin disk.

3

We neglect possible vertical variations in the gas temperature. A more accurate approach involving the reconstruction of the vertical structure is currently in development.

4

The adopted abundances are not expected to influence our results because we do not take chemical transformations into account in this study.

5

This value does not take the time spent at and Tmp > 20 K into account. The full time needed for the gas-phase fraction of CO to go back down to 0.5 following the burst is about 2−3 kyr.

Acknowledgments

The authors are thankful to Ruud Visser, the referee, for suggestions that helped to improve the manuscript and to Shu-ichiro Inutsuka for stimulating discussions. This work is supported by Royal Society awards WM090065 and RFBR Cost shared application with Russia (JP101297 and 11-02-92601). It was also partly supported by the European Research Council under the European Communitys Seventh Framework Programme (FP7/2007-2013 Grant Agreement No. 247060) and by the Consolidated STFC grant ST/J001627/1. The simulations were performed on the Shared Hierarchical Academic Research Computing Network (SHARCNET), on the Atlantic Computational Excellence Network (ACEnet), and on the Vienna Scientific Cluster (VSC-2).

References

  1. Armitage, P. J., Livio, M., & Pringle, J. E. 2001, MNRAS, 324, 705 [NASA ADS] [CrossRef] [Google Scholar]
  2. Baraffe, I., & Chabrier, G. 2010, A&A, 521, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  3. Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 2002, A&A, 382, 563 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  4. Baraffe, I., Chabrier, G., & Gallardo, J. 2009, ApJ, 702, L27 [NASA ADS] [CrossRef] [Google Scholar]
  5. Baraffe, I., Vorobyov, E. I., & Chabrier, G. 2012, ApJ, 756, 118 [NASA ADS] [CrossRef] [Google Scholar]
  6. Basu, S. 1997, ApJ, 485, 240 [CrossRef] [Google Scholar]
  7. Basu, S., & Mouschovias, T. Ch. 1994, ApJ, 432, 720 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 [NASA ADS] [CrossRef] [Google Scholar]
  9. Bisschop, S. E., Fraser, H. J., Oberg, K. I., van Dishoeck, E. F., & Schlemmer, S. 2006, A&A, 449, 1297 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Bonnell, I., & Bastien, P. 1992, ApJ, 401, L31 [Google Scholar]
  11. Bourke, T. L., Myers, P. C., Evans, N. J., II, et al. 2006, ApJ, 649, L37 [NASA ADS] [CrossRef] [Google Scholar]
  12. Chabrier, G., & Baraffe, I. 1997, A&A, 327, 1039 [NASA ADS] [Google Scholar]
  13. Chabrier, G., & Baraffe, I. 2000, ARA&A, 38, 337 [NASA ADS] [CrossRef] [Google Scholar]
  14. Charnley, S. B., Rodgers, S. D., & Ehrenfreund, P. 2001, A&A, 378, 1024 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2010, A&A, 510, L3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. D’Angelo, C. R., & Spruit, H. C. 2010, MNRAS, 406, 1208 [NASA ADS] [Google Scholar]
  17. Dapp, W. B., & Basu, S. 2009, MNRAS, 395, 1092 [NASA ADS] [CrossRef] [Google Scholar]
  18. Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [Google Scholar]
  19. Dunham, M. M., & Vorobyov, E. I. 2012, ApJ, 747, 52 [NASA ADS] [CrossRef] [Google Scholar]
  20. Dunham, M. M., Evans, N. J., II, Bourke, T. L., et al. 2006, ApJ, 65, 945 [NASA ADS] [CrossRef] [Google Scholar]
  21. Dunham, M. M., Evans, II, N. J., Terebey, S., Dullemond, C. P., & Young, C. H. 2010, ApJ, 710, 470 [NASA ADS] [CrossRef] [Google Scholar]
  22. Enoch, M. L., Evans, II, N. J., Sargent, A. I., & Glenn, J. 2009, ApJ, 692, 973 [NASA ADS] [CrossRef] [Google Scholar]
  23. Evans, N. J., II, Dunham, M. M., Jorgensen, J. K., et al. 2009, ApJSS, 181, 32 [Google Scholar]
  24. Harries, T. 2011, MNRAS, 411, 1500 [NASA ADS] [CrossRef] [Google Scholar]
  25. Harries, T., Monnier, J. D., Symington, N., & Kurosawa, R. 2004, MNRAS, 350, 565 [NASA ADS] [CrossRef] [Google Scholar]
  26. Hartmann, L., Zhu, Z., & Calvet, N. 2011 [arXiv:1106.3343] [Google Scholar]
  27. Hennebelle, P., & Teyssier, R. 2008, A&A, 477, 25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Ilee, J. D., Boley, A. C., Caselli, P., et al. 2011, MNRAS, 417, 2950 [NASA ADS] [CrossRef] [Google Scholar]
  29. Johnson, B. M., & Gammie, C. F. 2003, ApJ, 597, 131 [NASA ADS] [CrossRef] [Google Scholar]
  30. Kenyon, S. J., Hartmann, L. W., Strom, K. M., & Strom, S. E. 1990, ApJ, 99, 869 [Google Scholar]
  31. Kim, H. J., Evans, II, N. J., Dunham, M. M., et al. 2011, ApJ, 729, 84 [NASA ADS] [CrossRef] [Google Scholar]
  32. Kim, H. J., Evans, N. J., II, Dunham, M. M., Lee, J.-E., & Pontoppidan, K. M. 2012, ApJ, 758, 38 [NASA ADS] [CrossRef] [Google Scholar]
  33. Lee, J.-E. 2007, J. Korean Astron. Soc., 40, 83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Lin, D. N. C., & Papaloizou, J. 1986, ApJ, 309, 846 [NASA ADS] [CrossRef] [Google Scholar]
  35. Lodato, G., & Clarke, C. J. 2004, MNRAS, 353, 841 [NASA ADS] [CrossRef] [Google Scholar]
  36. Lucy, L. B. 1999, A&A, 344, L282 [NASA ADS] [Google Scholar]
  37. Masson, J., Teyssier, R., Mulet-Marquis, C., Hennebelle, P., & Chabrier, G. 2012, ApJS, 201, 24 [NASA ADS] [CrossRef] [Google Scholar]
  38. Machida, M. N., Inutsuka, S., & Matsumoto, T. 2011a, ApJ, 729, 42 [NASA ADS] [CrossRef] [Google Scholar]
  39. Machida, M., Inutsuka, S.-I., & Matsumoto, T. 2011b, PASJ, 65, 555 [NASA ADS] [CrossRef] [Google Scholar]
  40. Mathis, J. S., Rumpl, W., & Nordsieck, K. H. 1977, ApJ, 217, 425 [NASA ADS] [CrossRef] [Google Scholar]
  41. McKee, C., & Ostriker, E. C. 2007, ARA&A, 45, 565 [NASA ADS] [CrossRef] [Google Scholar]
  42. Noble, J. A., Congiu, E., Dulieu, F., & Fraser, H. J. 2012, MNRAS, 421, 768 [NASA ADS] [Google Scholar]
  43. Pfalzner, S., Tackenberg, J., & Steinhausen, M. 2008, A&A, 487, L45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Pontoppidan, K. M., Boogert, A. C. A., Fraser, H. J., et al. 2008, ApJ, 678, 1005 [NASA ADS] [CrossRef] [Google Scholar]
  45. Shu, F. S. 1977, ApJ, 214, 488 [NASA ADS] [CrossRef] [Google Scholar]
  46. Stamatellos, D., Whitworth, A. P., & Hubber, D. A. 2011, ApJ, 730, 32 [NASA ADS] [CrossRef] [Google Scholar]
  47. Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 134 [Google Scholar]
  48. Viallet, M., & Baraffe, I. 2012, A&A, 546, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Visser, R., & Bergin, E. A. 2012, ApJ, 754, 18 [NASA ADS] [CrossRef] [Google Scholar]
  50. Visser, R., van Dishoeck, E. F., Doty, S. D., & Dullemond, S. P. 2009, A&A, 495, 881 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  51. Vorobyov, E. I. 2012, Astron. Rep., 56, 179 [NASA ADS] [CrossRef] [Google Scholar]
  52. Vorobyov, E. I., & Basu, S. 2005, MNRAS, 360, 675 [NASA ADS] [CrossRef] [Google Scholar]
  53. Vorobyov, E. I., & Basu, S. 2006, ApJ, 650, 956 [NASA ADS] [CrossRef] [Google Scholar]
  54. Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896 [NASA ADS] [CrossRef] [Google Scholar]
  55. Zhu, Z., Hartmann, L., Gammie, C., et al. 2010, ApJ, 713, 1134 [NASA ADS] [CrossRef] [Google Scholar]
  56. Zhu, Z., Hartmann, L., Nelson, R. P., & Gammie, C. F. 2012, ApJ, 746, 110 [NASA ADS] [CrossRef] [Google Scholar]

All Figures

thumbnail Fig. 1

Gas surface density images in the inner 1000 × 1000 AU showing the disk evolution during 0.24 Myr after the formation of the central star. Vigorous disk fragmentation is evident in the figure. The scale bar is in log g cm-2.
In the text
thumbnail Fig. 2

Main characteristics of the forming system as a function of time elapsed since the formation of the central star: a) the mass accretion rate onto the protostar; b) the total and photospehric luminosities; c) the masses of the protostar, disk, and envelope; and d) the disk radius.
In the text
thumbnail Fig. 3

Effect of the luminosity bursts on the spatial distribution of gas-phase CO and CO2. The four columns present (from left to right) the gas surface density distribution (log g cm-2), gas temperature (K), CO gas-phase fraction, and CO2 gas-phase fraction in the inner 3000 × 3000 AU region. The four rows show (from top to bottom) the pre-burst phase (t = 0.128 Myr), the burst phase (t = 0.13 Myr), the post-burst phase (t = 0.134 Myr) and the quiescent phase (t = 0.14 Myr). The black lines are isotemperature contours outlining the regions where gas temperature exceeds 20 K (third column) and 40 K (fourth column) above which CO and CO2 are supposed to be in the gas phase in an equilibrium, time-independent case.
In the text
thumbnail Fig. 4

CO gas-phase fraction (red lines) and total stellar luminosity L∗ (black lines) vs. time. Two time periods (20 kyr each) are shown to highlight two strongest luminosity bursts and several moderate ones. The correlation between and L∗ is evident. In particular, steeply rises during the burst to a maximum value and gradually declines to a minimum value after the burst.
In the text
thumbnail Fig. 5

Top: CO gas-phase fraction vs. gas temperature Tmp calculated in the 500 < r < 1000 AU radial annulus during a time period covering a strong luminosity burst at t = 0.13 Myr. The dashed, solid, and dash-dotted lines show the pre-burst, burst, and post-burst phases with durations of 1.0 kyr, 0.2 kyr, and 5 kyr, respectively. The arrows indicate this evolution sequence. The horizontal and vertical dotted lines are the CO fraction of 0.5 and the gas temperature at which the CO desorbs. The grey-shaded area highlights the phase-space region with an abnormally high fraction of the gas-phase CO. The thick dash-dotted line highlights the –Tmp model track in this region. Bottom: CO2 gas-phase fraction vs. Tmp calculated in the 100 < r < 300 AU radial annulus during a strong luminosity burst at t = 0.082 Myr. The dashed, solid, and dash-dotted lines show the pre-burst, burst, and post-burst phases with durations of 0.4 kyr, 0.2 kyr, and 0.2 kyr, respectively. The horizontal and vertical dotted lines mark the CO2 fraction of 0.5 and the gas temperature at which the CO2 desorbs. The grey-shaded area highlights the phase-space region with an abnormally high fraction of the gas-phase CO2.

In the text
thumbnail Fig. 6

Top: gas radial temperature profiles. In particular, the red line presents the azimuthally averaged midplane temperature Tmp derived from hydrodynamical simulations at t = 0.131 Myr and the total stellar luminosity L∗ = 9  L⊙. The black solid line shows the gas temperature derived using the TORUS radiation transfer code, for the same gas density distribution and stellar luminosity as in the hydro model (case 1). The dashed solid line presents the TORUS calculation of the gas temperature for the same gas density distribution (case 1) but a higher stellar luminosity of 100  L⊙. Finally, the blue line shows the gas temperature calculated by the TORUS code for a spherically symmetric gas density distribution representing an infalling spherical envelope with the same mass as in the hydro model (case 2). Bottom: the gas volume density distribution used by the TORUS code (case 1) to calculate the gas temperatures (black solid and dashed lines). The scale bar is in g cm-3.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.