From Birth to Death of Protoplanetary Disks

From Birth to Death of Protoplanetary Disks: Modeling Their Formation, Evolution, and Dispersal


Formation, evolution, and dispersal processes of protoplanetary disks are investigated and the disk lifetime is estimated. Gravitational collapse of a pre-stellar core forms both a central star and a protoplanetary disk. The central star grows by accretion from the disk, and irradiation by the central star heats up the disk and generates thermal wind, which results in the disk dispersal. Using the one-dimensional diffusion equation, we numerically calculate the evolution of protoplanetary disks. To calculate the disk evolution from formation to dispersal, we add source and sink terms that represent gas accretion from pre-stellar cores and photoevaporation, respectively. We find that the disk lifetimes of typical pre-stellar cores are around 2–4 million years (Myr). A pre-stellar core with high angular momentum forms a larger disk whose lifetime is long, while a disk around a X-ray luminous star has a short lifetime. Integrating the disk lifetimes under various mass and angular velocity of pre-stellar cores and X-ray luminosities of young stellar objects, we obtain disk fraction at a given stellar age and mean lifetime of the disks. Our model indicates that the mean lifetime of protoplanetary disks is 3.7 Myr, which is consistent with the observational estimate from young stellar clusters. We also find that the dispersion of X-ray luminosity is needed to reproduce the observed disk fraction.

protoplanetary discs – accretion, accretion discs – circumstellar matter – stars: formation – stars: low-mass

1 Introduction

Protoplanetary disks are the birth place of planets. Revealing their formation, evolution, and dispersal processes is essential to understand planet formation.

The system of a protostar and a protoplanetary disk is formed by gravitational collapse of a rotating pre-stellar core in a molecular cloud. Observations provide us the statistical information of mass (e.g. André et al., 2014) and velocity gradient (e.g. Caselli et al., 2002) of pre-stellar cores. Using these information as an initial condition, three-dimensional numerical simulations of collapsing pre-stellar cores have unveiled the formation process of protostars and protoplanetary disks (e.g. Machida et al., 2010; Inutsuka, 2012; Tsukamoto et al., 2015). However, such three-dimensional simulations cannot follow the long-term evolution of the system of a protostar and a disk because of their expensive computational cost.

Most of gas in a pre-stellar core accretes onto a disk because of the centrifugal force. For the early evolutionary phase of this system, the disk is gravitationally unstable because the disk is more massive than the protostar (e.g. Vorobyov & Basu, 2006; Machida et al., 2010; Inutsuka et al., 2010). During this phase, gravitational torque by spiral arms formed in the disk efficiently transports the angular momentum, and the disk supplies the mass and entropy to the protostar through accretion, which affects evolutionary tracks of protostars (Baraffe et al., 2009; Hosokawa et al., 2011; Baraffe et al., 2012). After most of the gas accretes onto the disk, the protostar mass continues to increase whereas the disk mass decreases because of mass accretion from the disk to the protostar, so that the disk becomes gravitationally stable. Then, the turbulent stress caused by magneto-rotational instability (MRI) plays an important role for disk evolution (Balbus & Hawley, 1991). After the accretion timescale of the protostar becomes shorter than its thermal timescale, it starts the Kelvin-Helmholtz contraction (pre-main sequence, Hayashi, 1961).

The protoplanetary disks are observed as infrared excesses in the spectroscopic observation of young stellar objects (YSOs). Infrared observations of young star clusters give the relation between the cluster age and disk fraction. These data suggest that protoplanetary disks disappear in several million years (Myr) after the formation of pre-main sequence stars (e.g. Yasui et al., 2010; Mamajek, 2009). The widely accepted picture of disk dispersal is the combination of viscous accretion onto the star with the photoevaporation that is the thermally-driven wind induced by the irradiation of ultraviolet (UV) and X-rays (e.g. Clarke et al., 2001). Previous theoretical studies claim that the combination of photoevaporation and viscous evolution can explain the observed disk lifetime (e.g. Alexander & Armitage, 2009; Gorti et al., 2009; Owen et al., 2010; Bae et al., 2013). However, the formation process of disk has been ignored in these studies except for Bae et al. (2013).

The formation and dispersal processes of protoplanetary disks have been independently investigated so far. However, the disk evolution depends on the initial surface density profile, which is determined by the property of the pre-stellar core. Therefore, to obtain the reliable disk lifetime, the initial disk profile should be given by the model that takes account of collapse of a pre-stellar core. Constructing such a theoretical model enables us to bridge the observations of pre-stellar cores and disk lifetime. Bae et al. (2013) construct a model including the formation, evolution, and dispersal of the disk consistently. They calculate the evolution of the disk formed from pre-stellar cores of different angular velocities, and their model seems to be consistent with the disk lifetime estimated by observations. However, they do not consider the mass variation of pre-stellar cores, which is likely to affect the disk fractions estimated from observations of young star clusters.

In this paper, we calculate the formation, evolution, and dispersal processes of the protoplanetary disks. We start calculations from pre-stellar cores of various mass, and calculate the disk evolution with photoevaporation by X-rays until the disk disperses. We also estimate the disk fraction at a given stellar age and the mean lifetime of disks and compare it to the observations, which is the first attempt that bridges the observations of pre-stellar cores and disk lifetime. This paper is organized as follows. We describe a model of disk formation, evolution of a star and a disk, and photoevaporation in Section 2. The results of our calculations are represented in Section 3, and Section 4 is devoted to summary.

2 Model

We construct the model including disk formation from a pre-stellar core, co-evolution of a star-disk system, and photoevaporation by X-rays from the central star. To see a wide range of parameter space based on observations of pre-stellar cores, we use some simplifications compared to Bae et al. (2013).

2.1 Disk evolution

We solve the one-dimensional diffusion equation to calculate the evolution of protoplanetary disks. We use the cylindrical coordinate (), assume the axisymmetric flow, and integrate the fluid equations in vertical direction. Then, using the equations of mass and angular momentum conservation, we obtain the diffusion equation (e.g. Pringle, 1981),


where is the surface density, is the specific angular momentum, is the angular velocity, and is the kinetic viscosity. In addition to the viscous evolution, we consider the mass loss term by photoevaporation wind, , and the mass source term by infall from a pre-stellar core, , assuming that the specific angular momenta of the infall and wind are the same as those of the disk. We assume the hydrostatic equilibrium in the vertical direction, and then, the scale height is given as , where is the sound speed. We write , where is the Boltzmann constant, is the temperature of disk mid-plane, is the mean molecular weight, and is the mass of hydrogen atom.

In Equation (1), we assume the centrifugal balance, ( is the gravitational potential). The self-gravity of a disk should be included, since it is comparable to the stellar gravity in our calculations. We approximately write the gravitational potential as the spherically symmetric formula, , where is the gravitational constant and is the enclosed mass of the star and the disk inside the radius ( is the mass of the central star). Then, the specific angular momentum and the angular velocity are written as and , respectively.

We use the alpha prescription for viscosity (Shakura & Sunyaev, 1973). In protoplanetary disks, the gravitational torque and/or the torque due to magneto-rotational instability (MRI) are responsible for the angular momentum transport, which gives , where and represent the torque due to self-gravity and MRI turbulence, respectively. Since the gravitational torque depends on how the disk is gravitationally unstable, is often given as a function of Toomre’s parameter , where is epicyclic frequency (Toomre, 1964). We use according to Zhu et al. (2010) (see also Takahashi et al., 2013). For the turbulent viscosity by MRI, we use a constant value of alpha, . In reality, the strength of MRI turbulence depends on the resistivity, so that is expressed as a function of ionization degree (Gammie, 1996; Okuzumi & Hirose, 2011). Since this makes disk evolution complicated, we ignore the effect of resistivity as a first step study.

Since the thermal timescale is much shorter than the viscous timescale in the usual protoplanetary disks, we simply determine the temperature from the balance between three heating processes and a radiative cooling. As the heating processes, we consider viscous heating, accretion from the pre-stellar core, and irradiation from the central star and ambient molecular cloud. The viscous heating rate is


We assume that the specific internal energy brought by infalling material has the same value as the specific kinetic energy of the local Keplerian motion, in which the heating rate by infalling material from the pre-stellar core is written as


where we consider the rest of energy brought from the pre-stellar core is radiated away. The irradiation flux, , due to the central star of luminosity is estimated to be (e.g. Menou & Goodman, 2004)


where is the Stephan-Boltzmann constant, is the temperature determined only by the irradiation from the central star, and is the shielding factor. We fix (Hueso & Guillot, 2005) to avoid the numerical oscillation, such as thermal instability (see below). The irradiation heating rate is represented as (e.g., Menou & Goodman, 2004),


where is the temperature of ambient gas in molecular cloud and is the optical depth ( is the opacity).

The radiation cooling rate is represented as (Hubeny, 1990; Menou & Goodman, 2004)


where is the temperature at disk midplane. We calculate the equilibrium temperature of each process, and from , , and , respectively, and estimate the disk temperature to be


We note that in our model, the disk temperature is usually determined by one dominant process of the three heating processes. In this situation, error is small for the case with the sum of the three temperatures, rather than the average. In reality, depends on the temperature, but for simplicity, we use the optical depth calculated in the previous time step when calculating and ( is independent of ). Opacity in protoplanetary disk is given in Zhu et al. (2009) as a useful fitting formula. This opacity induces thermal instability for optically thick disks for K (Kimura & Tsuribe, 2012), which leads to the variable mass accretion onto protostar (Bell & Lin, 1994; Zhu et al., 2010; Ohtani et al., 2014). Since investigating the dynamics of thermal instability is beyond the scope of this study, we use the opacity of solid grains, , even if temperature is higher than the evaporation temperature of solid grains. This treatment cannot trace the variability of mass accretion rate in a short period (typically years). However, we consider that the time-averaged mass accretion rate is correctly calculated because the averaged mass accretion rate is determined at the outer region of the disk that is thermally stable. Owing to this treatment of thermal process, we can avoid shortening the timestep. This enables us to investigate a wide parameter range.

2.2 Infall from a pre-stellar core

To obtain , we use the model of disk formation based on Takahashi et al. (2013). We set a critical BE density profile with central density and sound speed as an initial state of the pre-stellar core. The radius of this core is , which is so called critical Bonnor-Ebert radius. Then, we increase the density as a factor of to make the gravity stronger than the pressure gradient. We assume that the core rigidly rotates with angular velocity .

Because of the density enhancement, the core starts gravitational collapse, and forms the disk and star system at the center of the pre-stellar core. We divide the pre-stellar core into spherical shells, and suppose that each shell of initial radius simultaneously accretes onto the disk and star. The time when the shell accretes is


where is the enclosed mass inside the shell of radius . The integral represents the delay of gravitational collapse due to the pressure gradient. The total mass accretion rate onto the disk and star is . We assume that the angular momentum of each fluid element conserves during the collapse, and the fluid element accretes at the radius where the gravitational force balances the centrifugal force. Then, the mass accretion rate at each radius is described as


The gas infall from the pre-stellar core to the disk ends at . For , and . See Takahashi et al. (2013, 2016) for more detail of the model of disk formation. This model reproduces a disk structure obtained by a three-dimensional hydrodynamical simulation of Machida et al. (2010).

As an initial condition of star and disk, we put a protostar and a disk of . We use the self-similar solution of Lynden-Bell & Pringle (1974) as an initial profile of surface density. Note that the initial profile is unimportant because the disk immediately forgets the initial condition due to the mass accretion from the pre-stellar core and viscous evolution. The time is the instant at which the core begins to collapse.

2.3 Photoevaporation

We consider photoevaporation by X-rays from the central star. Although the extreme ultra-violet (EUV) photons have also ability to generate disk winds, they cannot reach the dense region of the disk due to their large cross-section. Thus, the mass loss rate by EUV photons can be neglected compared to that by X-rays (e.g., Owen et al., 2010; Owen et al., 2012).

We use the total mass loss rate and mass loss profile provided by Owen et al. (2012). In their model, the total mass loss rate, , and mass loss profile, , are


where is the X-ray luminosity from the central star (see the next subsection for prescription of ) and is the normalization factor determined by the condition . We set for . In the late phase, this X-ray driven winds creates an inner hole in the disk. We evaluate the inner hole radius, , such that , where is the optical depth for X-rays from the central star ( cm is the X-ray cross-section for a particle of gas). The existence of an inner hole modifies the mass loss rate and profile. If is satisfied, we use the mass loss rate and profile with an inner hole written as


See Appendix B of Owen et al. (2012) for the detailed formulae of and .

After the hole formation, there remains a little amount of gas inside the hole. This gas can absorb the X-rays, so that the disk wind should blow even for . Since absorption rate of the X-rays is proportional to the density, we assume the following formula for ,


We set for .

It is likely that the photoevaporation is suppressed during the infall phase. Thus, for , we set .

Figure 1: The time evolution of physical quantities for run A2. The top panel shows (solid) and (dashed). The middle panel depicts mass accretion rate from disk to star (solid), (dashed), and (thick-dotted). The bottom panel represents (solid) and (dashed). The left and right thin-vertical dotted lines show and , respectively.

2.4 Evolution of central stars

We need to take the stellar evolution into account because its luminosity is important for the disk temperature. When the protostellar mass is increasing, we assume that the protostar evolves along the birthline provided by Stahler & Palla (2005), which corresponds to the evolutionary track of an accreting star. Although Baraffe et al. (2009) has pointed out that entropy brought to accreting stars through accretion can modify the evolutionary tracks of the accreting stars significantly, we use a standard evolutionary model for simplicity. When the growth of stellar mass is completed, we use the evolutionary tracks for pre-main sequence stars of Siess et al. (2000). We judge the end of growth of the stellar mass by the condition that the disk evaporation timescale becomes sufficiently shorter than the growth timescale of stellar mass, , where and ( is the mass accretion rate onto the central star). We define the instant when this condition is satisfied as the starting time of pre-main sequence phase, . We have confirmed that the growth of the stellar mass is less than after . We obtain the photospheric stellar luminosity, , and the stellar radius, , as a function of the stellar mass and time from the data sets of evolutionary tracks by using the linear interpolation. We also calculate the accretion luminosity, , and use the sum as the stellar luminosity .

We use the X-ray luminosity suggested by Preibisch et al. (2005): , where is a parameter in the present paper. The mean value and standard deviation of the parameter are set to be and dex, respectively.

2.5 Disk dispersal condition

The disk fraction of each cluster is usually determined by observations of near-infrared (NIR, e.g. Haisch et al., 2001), and the mean disk lifetime is derived by the integration of disk fractions of clusters. Thus, we use the disk dispersal condition such that the optical depth in the NIR-emitting region becomes lower than unity. We define the NIR-emitting region as a hotter region than 300 K because NIR ranges from 1 to 8 m. The disk dispersal time in NIR range, , is defined as the instant when is satisfied in the region where 300 K. Observationally, the age of young clusters is determined by using the evolutionary tracks of pre-main sequence stars. The age is not defined as the time after the collapse of a pre-stellar core, but defined as the time after the pre-main sequence phase starts. Thus, we estimate the disk lifetime of NIR to be .

3 Results

3.1 Typical evolution of a protoplanetary disk

[10 cm] [km s pc] []
A1 2.0 0.5 1.0
A2 (reference) 2.0 1.5 1.0
A3 2.0 2.5 1.0
B1 0.5 1.5 1.0
B2 8.0 1.5 1.0
C1 2.0 1.5 0.2
C2 2.0 1.5 5.0
Table 1: Parameter sets for each runs.

We calculate the evolution of protoplanetary disks, taking account of formation and evaporation. We set computational region from 0.01 to 10 au. We equally divide the region in logarithmic space with 100 meshes, and solve Equation (1) explicitly using the forward-time centered-space diferencing scheme (Press et al., 1992). We fix , cm s, K, and . We calculate several values of , , and as tabulated in Table 1.

Figure 2: The time evolution of radial profile of surface density (upper) and temperature (lower) for run A2.

Figure 1 shows the time evolution of physical quantities for run A2. There are three phases: infall phase , viscous evolution phase , and pre-main sequence phase . In the infall phase, the mass accretion rate from cloud core is as high as , and both the protostar mass and the disk mass increase. The gravitational torque dominates over the MRI turbulence for the mechanism of angular momentum transport during this phase. After infall finishes, the disk mass decreases while the protostar continues to grow up, owing to the high mass accretion rate from the disk onto the protostar. The MRI turbulence mainly drives this mass accretion, and the typical mass accretion rate is about , which is consistent with observation of T-Tauri stars. After a few Myr, the disk becomes sufficiently lighter than the central star, and the pre-main sequence phase starts because the condition is satisfied. For , the central star becomes fainter because of Kelvin-Helmholtz contraction. The mass accretion rate onto the central star continues to decrease, while the photoevaporation is efficient. Thus, the disk becomes lighter and lighter, and finally, the disk disperses.

The upper panel of Figure 2 shows the radial profile of surface density with some time snapshots for run A2. In the early phase of disk evolution, the disk spreads outward due to the diffusion nature of viscous accretion disks (e.g. Lynden-Bell & Pringle, 1974; Pringle, 1981). The cutoff radius (kind of disk radius; see Appendix for the definition of ) increases with time, and is larger than at as shown in Table 2, where is the centrifugal radius of outermost region of the pre-stellar core. Although the photoevaporation process has little effect on disk evolution during the infall and viscous evolution phases, it determines the disk structure in pre-main sequence phase. Since the mass loss profile by photoevaporation has a peak, the gap appears around 1 au (e.g. Owen et al., 2010; Bae et al., 2013). After the gap formation, the inner disk that is detached from the outer disk accretes onto the central star in its viscous timescale ( years), which is much shorter than other timescales, such as the disk evaporation timescale ( years). Thus, the gap becomes an inner hole. Then, the mass loss profile changes to that with an inner hole (see Section 2.3). This profile has a sharp peak at the inner edge of the disk. Therefore, the disk is dissipated from the inside to outside (Clarke et al., 2001; Owen et al., 2010).

The bottom panel of Figure 2 shows the radial profile of temperature for run A2. During the infall phase, viscous heating dominates over the other heating mechanisms inside the cutoff radius. The heating by the infalling material is always sub-dominant. After the infall finishes, the irradiation heating determines the temperature for outer region of the disk (e.g. au at Myr). As surface density decreases, the viscous heating rate is gradually decreasing, and the radius at which moves inward. After , the stellar evolution results in the decrease of the irradiation temperature with time. Note that the duration of the Hayashi phase of solar-mass stars is about 10 years.

Figure 3: The disk lifetime as a function of (upper), (middle), and (lower). The dotted lines show the dependence of disk lifetime on these parameters.
[Myr] [Myr] [Myr] [M] [M] [M] [10 au] [10 au] [10 au]
A1 0.16 0.73 1.6 0.90 0.73/0.82 0.17/0.071 6.5 0.14 1.9
A2 (reference) 0.16 2.2 2.5 0.90 0.53/0.76 0.37/0.12 6.5 1.3 6.8
A3 0.16 3.5 3.0 0.90 0.44/0.73 0.47/0.15 6.5 3.4 12
B1 0.32 5.1 3.3 1.8 0.80/1.3 1.0/0.44 13 8.9 15
B2 0.079 0.88 1.8 0.45 0.34/0.42 0.11/0.032 3.2 0.16 3.0
C1 0.16 5.3 10 0.90 0.53/0.82 0.37/0.074 6.5 1.3 15
C2 0.16 0.81 0.76 0.90 0.53/0.69 0.37/0.18 6.5 1.3 3.0
Table 2: Resultant quantities for each runs. For and , we write their values at both and . For , the values at is described.

We tabulate the resultant physical quantities in Table 2. All the calculations qualitatively have the same evolutionary characteristics described above. We find that the timescales are longer for lower , higher , and lower . The cores with higher have higher angular momenta. This makes the disk radius larger, which leads to the longer disk evolution timescales and . For lower , is longer, since is proportional to the free-fall time. The cores with lower have higher angular momenta owing to the larger initial disk radii, so that they have the longer and . The photoevaporation rate is higher for higher , which results in the shorter and . We also find that the disk lifetime has a power-law dependence on parameters, , , and as illustrated in Figure 3. Thus, the stellar X-ray luminosity has the strongest impact on the disk lifetime.

For runs with higher , lower , and lower , the duration of the accretion phase () seems to be longer ( 3–5 Myr) than expected (Evans et al., 2009; Maury et al., 2011). If this long accretion phase is real under a system with a high , the age spread of pre-main sequence stars in clusters (e.g. Hillenbrand, 2009) can be solved by the scatter of angular momenta of pre-stellar cores. It is also expected that the longer accretion phase arises from our treatment of stellar evolution, for which the pre-main sequence phase starts when the growth of the protostar stops. It is expected that the pre-main sequence phase starts when the timescale of Kelvin-Helmholtz contraction becomes shorter than the growth time of the central star. This condition is approximately written as , which is satisfied for (see the bottom panel of Fig. 1). In order to obtain a solid conclusion, we should simultaneously solve the evolution of the central star with time-dependent mass accretion, which is beyond the scope of this paper.

3.2 Disk fraction

For comparison of our results to observations, we derive disk fraction as a function of stellar age. To obtain the disk fraction, we consider the distribution functions for , and . From observations of pre-stellar cores, we can obtain the distribution function of and the core mass function (CMF; the distribution function of ). Caselli et al. (2002) observed velocity gradients of 26 pre-stellar cores. We use their data to obtain the distribution function of . We find that the distribution function is well-fitted by the normal function,


with km s ps and km s ps. We plot the observational data and fitting function of the distribution function in Figure 4, in which we see that the fitting formula matches the observation well. For CMF, we use CMF of Aquila given in Könyves et al. (2010), which is represented as a log-normal function,


where and are the peak and dispersion of the log-normal function, respectively. In addition to the distributions of pre-stellar cores, we consider the distribution of from X-ray observations of YSOs. Assuming the distribution function for is the log-normal, we can write it as


where erg s and (Preibisch et al., 2005).

Figure 4: The distribution function for . The thin line shows the observational data, and the thick line shows the fitting result.

Assuming that these parameters have no correlation, the disk fraction at the stellar age is represented as


where is the Heaviside’s step-function. We note that the age is the instant when the pre-main sequence phase starts. We perform calculations for a parameter range of , , and . We equally divide the grid of in linear space, while the grids of and are equally divided in logarithmic space. The comparison of resultant disk fraction to that estimated from observations is shown in the upper panel of Fig. 5. Our model gives the disk fraction that exponentially decreases with time , where 3.7 Myr is the mean lifetime of disks. This lifetime is consistent with the observational estimate of 3.8 Myr obtained by Yasui et al. (2014), although there still remains uncertainty in the observational estimate (the mean lifetime is 2.5 Myr according to Mamajek (2009)). Our model gives the slightly larger disk fraction at a given time. We discuss the reason in the next section.

To see the parameter that mainly broaden the disk lifetime, we calculate the disk fraction varying only one parameter and fixing the other parameters as their reference values (, , ). The result is shown in the lower panel of Fig. 5. Although the disk fractions are step-like because of the lack of resolution in the parameter space, we find that the X-ray luminosity mainly determines the broadening of disk fraction, whereas the angular velocity and mass of pre-stellar cores are sub-dominant. Therefore, further studies are highly encouraged to reveal the origin of the X-ray luminosities of YSOs (e.g. Flaccomio et al., 2003; Preibisch et al., 2005; Telleschi et al., 2007). This broadening of disk fraction by the dispersion in X-ray luminosities has been reported by Owen et al. (2011).

In this paper, we assume that the alpha parameter by MRI turbulence, , is constant. For the case with low , the ionization degree in the disk decreases, and MRI can be inactivated (Igea & Glassgold, 1999; Turner & Sano, 2008; Fujii et al., 2014). Since in the MRI inactive region is lower than that in MRI active region, it may be correlated with the X-ray luminosity. If this effect is included in our model, the model with weaker (stronger) has a lower (higher) value of , which results in longer (shorter) . This probably causes rapider decrease of the disk fraction in early phase and slower decrease in rate phase than that shown in Figure 5. The impact of dead zone on the mean lifetime, , remains as a future work.

Figure 5: Calculated disk fraction as a function of time. The upper panel shows comparison of our model to observational estimates. The our model (solid line) is roughly consistent with the observational estimate by Yasui et al. (2014) (dashed line). The dotted line is another observational estimate (Mamajek, 2009). The lower panel shows the effect of each parameter on broadening the lifetime of disks. The effects of rotation (thick-solid) and envelope mass (thick-dashed) is sub-dominant, while the X-ray luminosity (thick-dotted) mainly determines the disk lifetime. The model prediction shown in upper panel (thin-solid) and mean lifetime (thin-dotted) are also shown.

4 Summary

We have calculated the evolution of protoplanetary disks, taking account of formation from pre-stellar cores and dispersal by X-ray photoevaporation. A typical pre-stellar core collapses to form the star-disk system in 0.2 Myr. After infall from the pre-stellar core finishes, the star mass increases through accretion from the disk, while the disk mass gradually decreasing. Then, a few Myr later, mass accretion rate onto the star becomes sufficiently low, and pre-main sequence phase starts. It takes another 2–4 Myr for photoevaporation by X-rays to create a gap in disk at au. Soon after the gap formation, gas of the inner disk accretes to the star, and disk becomes unobservable in NIR. If a pre-stellar core has higher angular velocity or lower central density, its evolution timescale is longer than typical one. Also, the lower X-ray luminosity of the central star is, the longer the evolution timescale is. In addition, we have estimated the disk fraction at a given stellar age, constructing the distribution functions of pre-stellar cores and X-ray luminosities of YSOs by using the observational data. As a result, the mean lifetime of protoplanetary disks is Myr, which is consistent with the observational estimate. We have also found that the X-ray luminosity has the strongest impact on the disk lifetime in our model.

Finally, we discuss the effects of ignored processes, which may be the reason why our model slightly overestimates the disk fraction systematically. First, we ignore far ultraviolet (FUV) photons from the central star. Although the FUV photoevaporation rate in Gorti et al. (2009) is comparable to the X-ray photoevaporation rate in Owen et al. (2012), it is complex and still a work in progress (see Alexander et al., 2014, for a more detailed review). The large scale magnetic field also seems to affect the disk lifetime, which is ignored in our model for simplicity. The angular momentum in a collapsing pre-stellar core is transported outward by magnetic braking during the infall phase (e.g., Machida et al., 2011). This makes disk radius smaller, which results in shorter disk lifetime. Also, the large scale magnetic field drives disk winds via turbulent thermal pressure (e.g., Suzuki et al., 2010), or magnetic centrifugal force (e.g. Blandford & Payne, 1982). These winds also shorten the disk lifetime through the transport of angular momentum and mass. This effect of disk winds may be dominant process of disk dispersal if there are strong vertical magnetic fields, as discussed in Bai (2016). On the other hand, if we include the resistivity, some part of the disk becomes MRI inactive (dead zone, e.g. Gammie, 1996; Sano & Miyama, 1999, see Section 3.2). This probably makes disk lifetime longer, because the materials accumulate in the dead zone (Bae et al., 2013). In addition, as discussed in Section 3.1, our treatment of stellar evolution is incomplete. In reality, should be shorter, which probably makes the disk lifetime longer. More realistic calculations including these effects remain as a future work.


We thank Jaehan Bae, Chikako Yasui, and Shu-ichiro Inutsuka for useful comments. MK is supported by Grants-in-Aid from MEXT of Japan (Grant: 23244027).

Appendix A Estimation of cutoff radius

Consider a disk with the following surface density profile,


where is the cutoff radius. The enclosed disk mass inside is estimated to be


where is the lower incomplete gamma function. The total disk mass is described as


where is the gamma function. Dividing Equation (22) of by (23), we obtain the cutoff radius from following equation,


The disk profile is when the irradiation from the central star is the dominant heating source. In this case, . Using this value and (24), we tabulate the resultant cutoff radius at in Table 2. Independently, we get by directly fitting the simulation result by Equation (21), and it is found that these two estimates are consistent within a factor 1.2.


  1. pubyear: 2016
  2. pagerange: From Birth to Death of Protoplanetary Disks: Modeling Their Formation, Evolution, and DispersalA


  1. Alexander R. D., Armitage P. J., 2009, \apj, 704, 989
  2. Alexander R., Pascucci I., Andrews S., Armitage P., Cieza L., 2014, Protostars and Planets VI, pp 475–496
  3. André P., Di Francesco J., Ward-Thompson D., Inutsuka S.-I., Pudritz R. E., Pineda J. E., 2014, Protostars and Planets VI, pp 27–51
  4. Bae J., Hartmann L., Zhu Z., Gammie C., 2013, \apj, 774, 57
  5. Bai X.-N., 2016, \apj, 821, 80
  6. Balbus S. A., Hawley J. F., 1991, \apj, 376, 214
  7. Baraffe I., Chabrier G., Gallardo J., 2009, \apjl, 702, L27
  8. Baraffe I., Vorobyov E., Chabrier G., 2012, \apj, 756, 118
  9. Bell K. R., Lin D. N. C., 1994, \apj, 427, 987
  10. Blandford R. D., Payne D. G., 1982, \mnras, 199, 883
  11. Caselli P., Benson P. J., Myers P. C., Tafalla M., 2002, \apj, 572, 238
  12. Clarke C. J., Gendrin A., Sotomayor M., 2001, \mnras, 328, 485
  13. Evans II N. J., et al., 2009, \apjs, 181, 321
  14. Flaccomio E., Damiani F., Micela G., Sciortino S., Harnden Jr. F. R., Murray S. S., Wolk S. J., 2003, \apj, 582, 398
  15. Fujii Y. I., Okuzumi S., Tanigawa T., Inutsuka S.-i., 2014, \apj, 785, 101
  16. Gammie C. F., 1996, \apj, 457, 355
  17. Gorti U., Dullemond C. P., Hollenbach D., 2009, \apj, 705, 1237
  18. Haisch Jr. K. E., Lada E. A., Lada C. J., 2001, \apjl, 553, L153
  19. Hayashi C., 1961, \pasj, 13
  20. Hillenbrand L. A., 2009, in Mamajek E. E., Soderblom D. R., Wyse R. F. G., eds, IAU Symposium Vol. 258, The Ages of Stars. pp 81–94 (arXiv:0812.1261), doi:10.1017/S1743921309031731
  21. Hosokawa T., Offner S. S. R., Krumholz M. R., 2011, \apj, 738, 140
  22. Hubeny I., 1990, \apj, 351, 632
  23. Hueso R., Guillot T., 2005, \aap, 442, 703
  24. Igea J., Glassgold A. E., 1999, \apj, 518, 848
  25. Inutsuka S.-i., 2012, Progress of Theoretical and Experimental Physics, 2012, 01A307
  26. Inutsuka S.-i., Machida M. N., Matsumoto T., 2010, \apjl, 718, L58
  27. Kimura S. S., Tsuribe T., 2012, \pasj, 64
  28. Könyves V., et al., 2010, \aap, 518, L106
  29. Lynden-Bell D., Pringle J. E., 1974, \mnras, 168, 603
  30. Machida M. N., Inutsuka S.-i., Matsumoto T., 2010, \apj, 724, 1006
  31. Machida M. N., Inutsuka S.-I., Matsumoto T., 2011, \pasj, 63, 555
  32. Mamajek E. E., 2009, in Usuda T., Tamura M., Ishii M., eds, American Institute of Physics Conference Series Vol. 1158, American Institute of Physics Conference Series. pp 3–10 (arXiv:0906.5011), doi:10.1063/1.3215910
  33. Maury A. J., André P., Men’shchikov A., Könyves V., Bontemps S., 2011, \aap, 535, A77
  34. Menou K., Goodman J., 2004, \apj, 606, 520
  35. Ohtani T., Kimura S. S., Tsuribe T., Vorobyov E. I., 2014, \pasj, 66, 112
  36. Okuzumi S., Hirose S., 2011, \apj, 742, 65
  37. Owen J. E., Ercolano B., Clarke C. J., Alexander R. D., 2010, \mnras, 401, 1415
  38. Owen J. E., Ercolano B., Clarke C. J., 2011, \mnras, 412, 13
  39. Owen J. E., Clarke C. J., Ercolano B., 2012, \mnras, 422, 1880
  40. Preibisch T., et al., 2005, \apjs, 160, 401
  41. Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical recipes in FORTRAN. The art of scientific computing
  42. Pringle J. E., 1981, \araa, 19, 137
  43. Sano T., Miyama S. M., 1999, \apj, 515, 776
  44. Shakura N. I., Sunyaev R. A., 1973, \aap, 24, 337
  45. Siess L., Dufour E., Forestini M., 2000, \aap, 358, 593
  46. Stahler S. W., Palla F., 2005, The Formation of Stars
  47. Suzuki T. K., Muto T., Inutsuka S.-i., 2010, \apj, 718, 1289
  48. Takahashi S. Z., Inutsuka S.-i., Machida M. N., 2013, \apj, 770, 71
  49. Takahashi S. Z., Tomida K., Machida M. N., Inutsuka S.-i., 2016, preprint, (arXiv:1604.05432)
  50. Telleschi A., Güdel M., Briggs K. R., Audard M., Palla F., 2007, \aap, 468, 425
  51. Toomre A., 1964, \apj, 139, 1217
  52. Tsukamoto Y., Takahashi S. Z., Machida M. N., Inutsuka S., 2015, \mnras, 446, 1175
  53. Turner N. J., Sano T., 2008, \apjl, 679, L131
  54. Vorobyov E. I., Basu S., 2006, \apj, 650, 956
  55. Yasui C., Kobayashi N., Tokunaga A. T., Saito M., Tokoku C., 2010, \apjl, 723, L113
  56. Yasui C., Kobayashi N., Tokunaga A. T., Saito M., 2014, \mnras, 442, 2543
  57. Zhu Z., Hartmann L., Gammie C., 2009, \apj, 694, 1045
  58. Zhu Z., Hartmann L., Gammie C., 2010, \apj, 713, 1143
This is a comment super asjknd jkasnjk adsnkj
The feedback cannot be empty
Comments 0
The feedback cannot be empty
Add comment

You’re adding your first comment!
How to quickly get a good reply:
  • Offer a constructive comment on the author work.
  • Add helpful links to code implementation or project page.