Free Access
Issue
A&A
Volume 599, March 2017
Article Number A49
Number of page(s) 16
Section Stellar structure and evolution
DOI https://doi.org/10.1051/0004-6361/201628260
Published online 24 February 2017

© ESO, 2017

1. Introduction

Since the pioneering work of Hayashi (1961), the first phase of stellar evolution, the so-called pre-main-sequence or PMS, is generally considered with the following simple approach: a star is formed with a large radius and contracts isotropically, and a huge release of gravitational potential energy heats the interior and yields a similarly large luminosity, which ensures convection in most of the interior of the star. As the star heats up, thermonuclear reactions set in, contraction is slowed, and a radiative zone begins to grow. In standard models for the contraction of a solar mass star, the growth of the radiative zone starts after 2 million years and the star reaches the main sequence (with an outer convective zone that is only about 2.5% of the total mass of the present Sun) after about 40 million years. This long phase, in which the star possesses a very deep convective zone or is even fully convective, almost guarantees a homogeneous composition in the stellar interior owing to extremely fast mixing in convective zones. A large portion of the gravitational energy of the accreted material is supposed to be given to the star, ensuring a large radius and important luminosity (e.g., Stahler et al. 1980).

Contrary to that ideal picture, a number of studies have revealed that the PMS evolution can be strongly affected by the way material is accreted onto the star during the accretion phase. If material loses entropy before or during the accretion onto the star, it can grow from a small radius and avoid the large quasi-static contraction phase (e.g., Mercer-Smith et al. 1984; Palla & Stahler 1992; Hartmann et al. 1997; Baraffe et al. 2009; Hosokawa et al. 2011; Dunham & Vorobyov 2012). This has strong consequences for the stellar evolutionary tracks and therefore the inferred physical properties of the star (Baraffe et al. 2009, 2012; Hosokawa et al. 2011). Finally, it controls the level to which planet formation affects stellar surface composition, as suggested from observations (Ramírez et al. 2009, 2011; Meléndez et al. 2009) or theoretical models (Chambers 2010; Guillot et al. 2014).

In this article, we wish to understand what controls the evolutionary tracks on the PMS, explain apparent differences seen in published results, and obtain limits on physically plausible evolutionary tracks from observational constraints. In a subsequent article, we will attempt to understand what controls the evolution of radiative and convective zones in PMS stars.

This paper is organized as follows. In Sect. 2, we describe our physical model and computation method for simulating the PMS evolution including accretion. In Sect. 3, we examine how, and to what extent, energy losses during accretion control PMS evolution. We also identify deuterium as playing a leading role in this evolution phase despite its fast burning nature. In Sect. 4, we explore the dependence on the entropy of accreting matter. In Sect. 5, we compile the results of Sects. 3 and 4 to evaluate observational constraints and observational consequences. Our results are summarized in Sect. 6.

2. Method

2.1. Stellar evolution code

Our evolution models are calculated for spherically symmetric single stars without strong magnetic fields or rotation, but including accretion from a small initial seed. We use the one-dimensional stellar evolution code MESA version 6596 (Paxton et al. 2011, 2013, 2015) and refer to the Paxton et al. papers for full details of the computational method.

The code numerically solves the equations of continuity, momentum, energy, temperature gradient, and composition. We assume the hydrostatic equilibrium for the momentum equation. The temperature gradient is determined by the mixing length theory of Cox & Giuli (1968) and Henyey et al. (1965). The Ledoux criterion of convection is used. The composition is changed by thermonuclear reaction and diffusion. The diffusion coefficient is given by the mixing length theory and overshooting of Herwig (2000), in which the diffusion coefficient exponentially approaches zero in the radiative zone from the boundary with the convective zone. We use the mixing length parameter, αMLT = 1.90506, which is the ratio of the mixing length to the pressure scale height, and the overshooting parameter, fov = 0.0119197, which determines the extent of the overshooting region (see Eq. (9) of Paxton et al. 2013) (see Appendix A). The energy equation is described in Sect. 2.6 in detail.

We use the radiative opacity of Seaton (2005, OP) and Ferguson et al. (2005) for the temperature range below 104.5 K, and the electron conduction opacity of Cassisi et al. (2007). We adopt the MESA default version of the equation of state; we basically follow the OPAL EOS tables (Rogers & Nayfonov 2002) and SCVH tables (Saumon et al. 1995).

2.2. Thermonuclear reactions

Pre-main-sequence stars ignite light elements such as deuterium and lithium when the central temperature exceeds ~ 106 and ~ 3 × 106 K, respectively. The deuterium burning is as follows: D + 1H → 3He + γ. This is a strongly exothermic reaction, which produces 5.494 MeV via one reaction. Although the lithium burning produces more energy by a factor of about three, it does not affect the evolution owing to a ~ 103 times smaller abundance (Grevesse & Sauval 1998). We use the thermonuclear reaction rates of Angulo et al. (1999, NACRE). Main-sequence (MS) stars ignite hydrogen when temperature exceeds ~ 107 K.

2.3. Chemical composition

The initial mass fractions of hydrogen, helium, and metals are denoted by Xini,Yini, and Zini, respectively. We choose the input parameters that reproduce values of the present-day Sun estimated by helioseismic and spectroscopic observations (see Appendix A). We conducted a χ2 test and used the following converged input parameters: Xini = 0.70046, Yini = 0.27948 and Zini = 0.02006. We assume 3He /4He = 10-4, similar to the value in the Jovian atmosphere (Mahaffy et al. 1998). We use the composition of metals of Grevesse & Sauval (1998).

As described in the previous section, the deuterium content is a key parameter for the PMS evolution. In this study we use the mass fraction of deuterium, XD, and our fiducial XD is 20 ppm (parts per million, 10-6), which is an interstellar value. However, the interstellar XD remains uncertain. Thus we explore the consequences of varying the deuterium content in Sect. 3.3.

Finally, the accreting gas is assumed to have constant composition with time. Importantly, this includes fresh deuterium, which is added to the outer layers. Since deuterium is quickly burned in the stellar interior, this is an important factor governing the evolution on the PMS. The deuterium abundance profile is calculated by accounting for diffusion.

2.4. Initial conditions

Stars are formed by the collapse of molecular cloud cores through a number of stages (e.g., Larson 1969; Winkler & Newman 1980; Inutsuka 2012). As the collapse of a molecular cloud core proceeds, the central density increases and the central region becomes adiabatic. After the first hydrostatic core is formed, the dissociation of hydrogen molecules causes the second collapse after the central temperature exceeds ~ 2000 K. As a result of the second collapse, the central temperature increases and eventually a second hydrostatic core is formed.

In this paper, we focus on the evolution after the formation of the second core. This second core is only about a Jovian mass (MJup) initially, and most of the mass is thus accreted in a subsequent phase. As accretion is progressively suppressed, the star enters its quasi-contraction phase, increases its central temperature, and eventually becomes a main-sequence star.

In this paper, we set our fiducial values of the initial mass and radius as 0.01 M and 1.5 R, respectively, following Stahler et al. (1980) and Hosokawa et al. (2011, hereafter HOK11. This choice is essentially driven by convergence issues with MESA at low masses. However, our fiducial initial mass is higher than the second core mass, which is derived by radiation-hydrodynamic simulations (~ 14 MJup; e.g., Masunaga & Inutsuka 2000; Vaytet et al. 2013). Therefore initial seeds in the present work correspond to slightly evolved protostars from second cores. As discussed by Baraffe et al. (2012, hereafter BVC12, the values of the initial seed mass and radius set the entropy of the forming star and therefore its subsequent evolution. We show in Appendix D that for ranges of masses between 1 and 4 MJup, the range of radii when the mass attains 0.01 M is between 0.25 R and 1.5 R. Our fiducial value of 1.5 R therefore corresponds to a high initial value of the entropy of the star. In Sect. 3.4 we explore the consequences of variations of the initial radius down to 0.25 R.

2.5. Mass accretion

The collapse of the molecular cloud core yields an accretion rate ~cs3/G\hbox{$\dot{M} \sim c\sub{s}^3/G$}, where cs is the characteristic speed of sound in the cloud (e.g., Shu 1977; Stahler et al. 1980). For molecular cloud temperatures T ~ 30 K, this implies ≈ 10-5M/ yr, which is our fiducial accretion rate. The accretion onto the star is however mostly determined by the angular momentum of the collapsing gas and formation of a circumstellar disk (e.g., Hueso & Guillot 2005; Vorobyov & Basu 2010; Inutsuka et al. 2010; Martin et al. 2012). It can be strongly variable (episodic), as observed in FU Ori.

Our fiducial accretion rate is ≈ 10-5M/ yr but we explore the effect of varying this rate in Sect. 3.2. We stop mass accretion abruptly when the stellar mass reaches the final mass, Mfin.

2.6. Modeling the consequences of accretion

The consequence of accretion onto the forming star is of course an increase of its mass, but this accretion also modifies the stellar environment and its radiation to space and delivers energy to the star. Accretion also delivers fresh deuterium, which is an important source of combustible material.

It is difficult at this point to model the problem in its full complexity, in particular, because this should involve accounting for the angular momentum of material in the molecular cloud core, its magnetic field, the presence of outflows, and the particular geometry that arises from the birth of a circumstellar disk. Instead, we adopt a simple parametric approach inspired by the work of Baraffe et al. (2009, hereafter BCG09.

First, we assume that the specific entropy of the accreting material is the same as that at the stellar surface. This follows from Hartmann et al. (1997) estimate that, in the case of accretion from disk to protostar, any entropy (equivalent temperature) excess would be quickly radiated to space because, especially in the presence of an inner cavity, the optical thickness toward the stellar interior is much larger than toward the disk. In practice, this treatment is favorable to the computational convergence. We point out that validating this hypothesis would require a multidimensional radiation-hydrodynamic simulation beyond the scope of the present work.

Second, we parametrize the heat injected by the accreting material as Ladd=ξGM/R,\begin{equation} L\sub{add}=\xi GM_\star\dot{M}/R_\star, \label{eq:Ladd} \end{equation}(1)where M and R are stellar mass and radius, respectively, and energy constraints impose that ξ is such that 0 ≤ ξ ≤ 11. In situations in which radiation to space is favored (e.g., when accreting from a disk with an inner cavity), we expect a low value of ξ. Even when this is not the case, Stahler et al. (1980) show that radiative transfer considerations impose that ξ ≤ 3/4 in the spherical accretion. If materials accrete onto the star through an active Keplerian disk, the upper limit of ξ should be 0.5 owing to the radiative cooling from the disk surface. If the star rotates rapidly, a large fraction of gravitational energy can be stored in the rotational energy and then ξ is decreased. Therefore, in practice we opt for ξ = 0.5 as a realistic upper limit for what we call “hot accretion”. We label simulations with ξ = 0 as “cold accretion” simulations and we label as “warm accretion” simulations all cases where 0 <ξ< 0.5.

Table 1

Input parameters.

Third, we use two simple models to parametrize how the energy enters the stellar interior. The first is the uniform model, which is adopted by Baraffe & Chabrier (2010, hereafter BC10, and in which Ladd is distributed uniformly and instantaneously within the entire star. The energy deposited per unit mass is thusεadd(uniform)=LaddM·\begin{equation} \varepsilon\sub{add}^{\rm (uniform)}=\frac{L\sub{add}}{M_\star}\cdot \label{eq:eps_uniform} \end{equation}(2)The second is the linear model, in which we consider that the accretion energy is deposited only in an outer layer of relative mass mke, such that the accretion energy is zero in the layer from the stellar center to relative mass (1−mke) and that it increases linearly with mass until the photosphere. With this assumption, the energy deposited per unit mass is εadd(linear)=LaddMMax[0,2mke2(MrM1+mke)]\begin{equation} \varepsilon\sub{add}^{\rm (linear)}=\frac{L\sub{add}}{M_\star}{\rm Max}\left[0,\frac{2}{m\sub{ke}^2}\left( \frac{M_{r}}{M_\star}-1+m\sub{ke} \right) \right] \label{eq:eps_linear} \end{equation}(3)where Mr is the mass coordinate and mke is expressed as a fraction of the total mass of the star at time t so that 0 <mke ≤ 1. In both cases, εadd satisfy the relation εadddM=Ladd\hbox{$\int \varepsilon\sub{add} {\rm d}M=L\sub{add}$}. The energy conservation equation thus is written ∂LMr=εnuc(T∂S∂t)Mr+εadd,\begin{equation} \frac{\partial L}{\partial M_{r}} = \varepsilon_{\mathrm{nuc}}-\left( T\frac{\partial S}{\partial t} \right)_{M_r}+\varepsilon_{\mathrm{add}}, \label{eq:dLdM} \end{equation}(4)where L is the luminosity, T its temperature, S its specific entropy, t time, and εnuc is the energy production rate from thermonuclear reactions in the shell.

The uniform model thus corresponds to the case studied by BC10 in which Ladd is distributed uniformly and instantaneously within the entire star. In the linear model, we consider that the accretion energy is deposited preferentially in the outer layers of the star and that this energy deposition is linear in mass over a shell of mass mke (with mass expressed as a fraction of the total mass of the star at time t).

BC10 justify the uniform model by the fact that low-mass (≲ 2 M) stars are fully convective during their PMS stage and that convection may rapidly deliver the accretion energy into the deep interior. We believe that its efficient transport of energy (or equivalently entropy) is probably unrealistic owing to the difficulty in transporting energy into a star. Furthermore, the generation of heat at the surface of accretion flows can potentially hinder convection, thereby preventing any uniform entropy mixing, as recently found by two-dimensional hydrodynamic simulations with high accretion rates (Geroux et al. 2016). The linear model is hence a highly simplified but useful parametrization that allows us to explore the consequences of a penetration of the accretion energy only to a fraction of the radius of a star.

In the case of hot accretion, HOK11 adopt a different approach; they do not consider that accretion energy can be transported inward but instead adopt an outer boundary condition to account for the heat generated by accretion (see Hosokawa & Omukai 2009). In principle, this should be a better approach. However, the approach remains one-dimensional and is not parametrized. It does not include the possibility that part of the energy may be transported to deeper levels in an non-radiative way. We therefore believe that the different approaches are complementary.

Numerically, evaluating the thermodynamic state of the accreting material requires careful consideration. The mass increase is performed within MESA using an “Eulerian” scheme (Sugimoto & Nomoto 1975). MESA originally treated the thermodynamic state of the accreting material with compressional heating (Townsley & Bildsten 2004; Paxton et al. 2013). However, this is not suitable for rapid accretion whose timescale is shorter than the thermal relaxation timescale. Thus, the new scheme was implemented in version 5527 and is used in this study. The detailed treatment is described in Sect. 7 of Paxton et al. (2015).

2.7. Outer boundary condition

The pressure Ps and temperature Ts at the outer boundary is specified by the interpolation of a model atmosphere which is calculated with the assumption that materials accrete onto a small fraction of the stellar surface from the thin disk and do not affect the properties of the photosphere.

We found that the choice of the model atmosphere table directly affects the convergence of the calculations. We hence selected the following two tables depending on stellar mass and radius:

For stars such that M ≥ 1 M and R> 0.7 R, we adopted a “photospheric” table for an optical depth τs = 2/3 and Ts = Teff and a surface pressure Ps set by the model atmospheres PHOENIX (Hauschildt et al. 1999a,b) and ATLAS9 (Castelli & Kurucz 2004). For less massive or smaller stars, we used the “τs = 100” model, which specifies Ps and Ts at τs = 100 from the ATLAS9 and COND (Allard et al. 2001) model atmospheres. The impact of the switch of the boundary condition on the stellar structure and evolution is negligible.

2.8. Adopted parameters and comparison with previous studies

Table 1 summarizes the values of our main parameters. Our fiducial values are shown in boldface. We also provide the range of values that we consider in Sects. 35 when we seek to explore the consequences of deviations of some parameters from the fiducial values. In particular, we explore the dependence on the deuterium mass mixing ratio, something which had not been considered so far.

We compare our results with two particular evolutionary models from previous studies (BC10 and HOK11) in Sect. 3.6.

3. Evolution in the context of a “cold” accretion

In this section we describe the resulting PMS evolution in the context of a “cold” accretion of material, i.e., in the extreme case when all of the accretion energy is radiated away (i.e., ξ = 0 in Eq. (1)).

Figure 1 shows how the radius changes with time in the cold accretion case and with our fiducial settings (see Table 1). As described in more detail in Appendix B, the evolution can be split into five phases: (I) the contraction phase; (II) the deuterium-burning phase; (III) the second contraction phase; (IV) the swelling phase; and (V) the main sequence. The contraction seen in phase (I) is induced by mass accretion with the same entropy as the star. Phase (II) begins when the internal temperature becomes high enough to start burning pre-existent deuterium, resulting in an expansion. With a decrease in the deuterium content, we enter phase (III) in which mass accretion again dominates over deuterium burning and the star contracts. After mass accretion is completed, the star expands owing to a change of its internal entropy profile. Eventually, after ~ 2×107 yr, it enters the main sequence. Starting from a different initial condition (e.g., a smaller initial radius) does not change the sequence, but it alters the duration and characteristics of the star before it reaches the main sequence.

After presenting the main differences between classical evolution models and those in the context of cold accretion (Sect. 3.1), we examine the influence of mass accretion rate, deuterium content, final stellar mass, and initial conditions in Sects. 3.23.5. We compare our results against previous studies in Sect. 3.6.

thumbnail Fig. 1

Evolution of radius (solid line) and mass (dashed) in the cold accretion limit (ξ = 0) and under the fiducial conditions listed in Table 1. The five main evolution phases are labeled as follows: (I) the contraction phase; (II) deuterium-burning phase; (III) second contraction phase until the accretion is completed; (IV) swelling phase; and (V) main sequence from ~ 2 × 107 yr.

3.1. Cold versus classical models

We now compare the radius evolution in the case of cold steady accretion with that in the classical model. Figure 2a shows that the stellar radius in the case of low-entropy accretion is one order of magnitude smaller than in the classical evolution case. Classically, a spherically symmetrical accretion prevents radiative losses implying a large stellar entropy and therefore a large radius. Following Stahler & Palla (2005), we adopt for the classical evolution case an initial radius of 4.92 R at 105 yr. When radiative cooling of accreting materials is allowed, the resulting specific entropy is much smaller, which yields a smaller protostar. Moreover, in the case of the low-entropy accretion, the MS starts at about 20 million years, which is faster than the classical evolution by a factor of about two. This difference comes from the smaller thermal timescale in the stellar interior because of the smaller radius in the case of low-entropy accretion.

thumbnail Fig. 2

Radius evolution (top panel) and evolutionary tracks in the H-R diagram (bottom) comparing the cold accretion (solid line) with the classical 1 M evolution (dashed). The filled square indicates the initial location, the triangle indicates the stellar mass at 0.1 M, and the circle indicates the end of accretion. The non-accreting classical evolution, whose initial radius is assumed to be 4.92 R as shown by the asterisks, starts at 0.99 × 105 yr to compare the evolution. At 105 yr, the radii are about one order of magnitude different. In the classical model the star evolves along the Hayashi and Henyey tracks, whereas in the cold accretion model the evolution is roughly horizontal in the H-R diagram. Although the deuterium burning slightly increases the luminosity where log Teff ≃ 3.5, it is much smaller than the Henyey track of the non-accreting 1 M star.

The pre-main-sequence evolutionary tracks for a 1 M star in the classical, non-accreting case and in the cold accretion scenario are shown in Fig. 2b. The tracks are drastically different. In the classical case, the star starts with a high entropy, then evolves along the almost vertical Hayashi track until it reaches the horizontal Henyey track in the Herzsprung-Russell (H-R) diagram. In the cold accretion scenario, the initial evolution in the H-R diagram is nearly horizontal because the accretion does not inject energy and both the entropy and intrinsic luminosity remain small. The implications of these results, including the dependence on ξ, are discussed in Sect. 5.

We point out that the evolutionary track in the cold accretion case in Fig. 2b does not include the accretion luminosity. Once the star moves out of the class I phase and becomes a visible T-Tauri star, its accretion luminosity is easily distinguished from its intrinsic luminosity in the spectral energy distribution because it is emitted in the UV or X-rays rather than in the visible. However, in the embedded phase, the accretion luminosity is reabsorbed by the surrounding molecular cloud and re-emitted at longer wavelengths, making the distinction with the intrinsic luminosity impossible. In that case, its location in the H-R diagram should be shifted upward (HOK11).

In the following three sections, we explore the consequences of changing the accretion history, deuterium content, initial entropy, and final mass. The other parameters are set to be the fiducial values listed in Table 1.

3.2. Dependence on mass accretion rate

As described in Sect. 2.5, the accretion rate onto the star can be highly time dependent and variable from one star to the next. We calculate the evolution varying mass accretion rate, ranging from 10-4 to 10-6M/ yr, and time variability; i.e., the episodic accretion following the pioneering works of BCG09, BC10 and HOK11. For episodic accretion, we adopt the parameters of BC10 on the basis of hydrodynamic simulations by Vorobyov & Basu (2005).

We confirm the findings of BCG09 and HOK11 that the variation of the accretion rate only has a moderate impact on the evolution. The difference in radii for the same mass is at most 0.2 R even for accretion rates that differ by two orders of magnitude. As described afterward, this is much smaller than the difference from other effects (the deuterium content XD and the heat injection efficiency ξ). The evolutionary tracks are also hardly affected by the variation of the accretion history. We stress that rather than accretion history, it is the entropy of the accreted material that matters. We therefore choose to fix the mass accretion rate to 10-5M/ yr from now on.

3.3. Dependence on deuterium content

The PMS evolution is controlled by deuterium burning when it occurs vigorously (i.e., during phase II in Sect. 3). The deuterium content can largely differ from star to star because of galactic evolution or the local environment. We now explore how PMS evolution is affected by deuterium content.

This was already investigated by Stahler (1988) who concluded that the radius is only moderately affected by deuterium content. For example, when stars become 1 M the radii are only increased by a factor of ≃ 1.6 from a deuterium-free star to the case with XD = 35 ppm. However, since they assumed a spherical accretion, this study is to be re-examined in the framework of low-entropy accretion.

First we summarize the currently available values of the deuterium content2. The primordial value at the beginning of the universe is XD,prim = 35.8 ± 2.5 ppm (Steigman 2006), the indirectly3 estimated value of the protosolar nebula is XD,PSN = 28.0 ± 2.8 ppm (Asplund et al. 2009), and the present-day value of the local interstellar medium varies widely from XD,ISM = 13.7 ± 5.3 ppm (Hébrard et al. 2005) to at least 32.3 ± 3.4 (Linsky et al. 2006). Classical studies, such as Stahler et al. (1980), used a higher value (XD = 35 ppm) based on classical observations of local interstellar media (e.g., Vidal-Madjar & Gry 1984).

Although these values still remain uncertain to some extent and, in particular, XD,ISM is still under debate (see Steigman 2006; Linsky et al. 2006; Prantzos 2007), they suggest that the deuterium content evolves with time. Considering that XD,PSN = 28 ppm at 4.57 Gyr ago, we assume that the deuterium content of present-day star-forming regions may be as low as 20 ppm and use that as our fiducial value. In addition to the time evolution, the deuterium content may show a wide range of values depending on the environment. For example, for stars formed from a cloud strongly affected by winds from nearby evolved stars, we could expect a smaller deuterium content than for other stars. It is therefore important to examine how variations in the deuterium fraction affect the PMS evolution of stars. In order to do so, we use the fiducial settings but vary the D/H ratio.

Figure 3 shows the radius evolution and evolutionary tracks with the different deuterium mass fraction from zero to 40 ppm, where the upper limit is set by the cosmic primordial value. Whereas the evolution before deuterium burning is totally independent of the deuterium content, we see that the tracks deviate after deuterium fusion sets in. For example, in the case XD = 40 ppm, the radius and luminosity are up to ~ 4 and ~ 30 times larger, respectively, than for the deuterium-free case. With the largest deuterium content, Lnuc exceeds 17Lacc\hbox{$\frac{1}{7}L\sub{acc}$} for a longer period of time and then the star becomes larger (see Eq. (B.4)), where LaccGM/R is the gravitational energy of the accreting material and Lnuc is the energy production rate by thermonuclear reaction. Evolutionary tracks are also largely different depending on XD in the region where log Teff ≳ 3.48, i.e., Teff ≳ 3000 K. These large differences illustrate the importance of the deuterium content on the PMS evolution.

thumbnail Fig. 3

Radius evolution (top panel) and evolutionary tracks (bottom) with varying deuterium contents ranging from XD = 40 (red solid line), 30 (green dashed), 20 (blue dotted), 10 (pink dot-dashed), and 0 (blue double dot-dashed) ppm. The points represent the end of the accretion phase and the black solid lines represent the classical PMS evolutionary tracks.

In the extreme case XD = 0, the star keeps shrinking even after ~ 104 yr. Just before the accretion is completed at 105 yr, the radius becomes ~ 0.3 R and eventually the star slightly expands. This expansion results from hydrogen burning due to sufficiently high central temperature even before the accretion is completed.

We thus find that the deuterium content affects low-entropy accretion evolution more strongly than found by Stahler (1988) with spherical accretion. This is because the total energy injected by accretion, Eacc=Ladddt\hbox{$E\sub{acc}=\int L\sub{add} \mathrm{d}t$}, exceeds the total energy generated by deuterium burning, ED = qDXDMfin, where qD = 2.63 × 1018 erg g-1 is the energy released by deuterium fusion per gram of deuterium, when ξ0.05(XD20ppm)(R)(MfinM)-1,\begin{equation} \xi \ga 0.05 \brafracket{X\sub{D}}{20~\mathrm{ppm}}\brafracket{\bar{R}}{\Rsun}{\brafracket{M\sub{fin}}{\Msun}^{-1}}, \end{equation}(5)where \hbox{$\bar{R}$} is the time-averaged stellar radius. In spherical accretion, Eacc always dominates ED and the difference in deuterium content is less pronounced. The reverse is true in the case of cold accretion, implying changes in the radii by up to a factor 3 and in luminosity by up to 2 orders of magnitude between the extreme XD values.

3.4. Dependence on initial conditions

So far we have assumed that our initial stellar seed is characterized by 1.5 R and 0.01 M. As pointed out by HOK11, since the entropy of accreting materials is related to the stellar entropy in our prescription (see Sect. 2.6), the initial conditions affect the evolution (see also Hartmann et al. 1997). This is also pointed out by BVC12 who stress the importance of the initial seed mass in the subsequent PMS evolution. Here, for computational reasons, we choose to vary the initial entropy (i.e., radius) instead of the initial mass. As discussed in Appendix D, the approach is equivalent with the current uncertainties in stellar seed masses leading to values of the radius at 0.01 M between 0.25 R and 1.5 R.

In Fig. 4, we show the evolutionary tracks of accreting protostars with brown-dwarf masses obtained with Rini from 0.25 R to 3 R and XD = 20 ppm and 35 ppm. We do not show the results for Rini ≳ 3 R since these are similar to the Rini = 3 R case because of their short K-H timescale. Like HOK11 (see their Figs. 5–7), we confirm that different Rini values can produce a scatter in the low-temperature region (by up to a 1.5 dex difference in luminosity), which cannot be produced by a varying deuterium abundance. Separately, we note that changing the initial conditions results in deuterium fusion setting in at a different stellar mass (0.08 and 0.03 M for Rini = 1.5 and 0.25 R, respectively; see Fig. D.1).

thumbnail Fig. 4

Evolutionary tracks for different initial radius Rini and different initial deuterium content XD assuming cold accretion (ξ = 0). The solid, dashed, dotted and dot-dashed lines represent the cases Rini = 0.25,0.5,1.5, and 3 R, respectively. The colors indicate deuterium content, XD = 20 ppm (blue) and XD = 35 ppm (red). The points represent the observed PMS stars in clusters (see Fig. 8). In all cases, stellar masses reach about 0.5 M at 3500 K.

We see that the evolutionary track with Rini = 0.25 R and XD = 20 ppm is characterized by a very low luminosity. The main reason for the small luminosity is the small entropy of both an initial seed and accreting materials. In the cold accretion cases, the entropy of the accreting material is assumed to be the same as the stellar entropy (see Sect. 2.6). Moreover, the effect of deuterium burning has less impact in this case. The star expands in Teff ≃ 25002900 K because of deuterium burning and luminosity increases. However the duration of the expansion phase is shorter than the other cases because of the large accreting material’s gravitational energy, Lacc4. As described in Appendix B, stars expand when LD>17Lacc\hbox{$L\sub{D}>\frac{1}{7}L\sub{acc}$} (Eq. (B.4)). Hence the larger Lacc makes the expansion phase (phase II in Sect. 3) shorter and the stellar radius and luminosity remain small (see also BVC12).

3.5. Evolution tracks for different final masses

In this section we focus on the evolution of stars with varying masses. Apart from the final masses, the settings are the same as the fiducial values listed in Table 1. The cold accretion evolution of stars with masses of 0.05, 0.1, 0.3, 1, and 1.5 M are shown in Fig. 5. Low-mass stars contract for a long time because of their long K-H timescale. For example, stars with 0.1 and 0.3 M enter their MS at 109 and 3 × 108 yr, respectively. During the first million years of evolution, Fig. 5 shows that the size of stars undergoing cold accretion is not a monotonic function of mass. For example, a 0.1 M star expands just after accretion ceases due to deuterium burning and becomes larger than the other stars considered here for some time. On the other hand, a 0.05 M star contracts monotonically and even more after 1–2 million years when all deuterium has been consumed. Higher mass stars (1.0, 1.5 M) expand just after accretion has stopped because of hydrogen burning.

thumbnail Fig. 5

Radius evolution (top panel) and evolutionary tracks (bottom) with varying final masses ranging from 0.05 M (solid line), 0.1 M (dashed), 0.3 M (dotted), 1 M (dot-dashed), and 1.5 M (double dot-dashed). The accretion is cold (ξ = 0) and steady (10-5M/ yr). The filled circles indicate the ages when the accretion is completed and the pluses, crosses, and asterisks represent 1, 5, and 10 million years, respectively. The calculation stops at 1010 yr or when stars start leaving the MS (only in the case of 1.5 M). We note that the evolutionary track of 0.05 M is cut short.

3.6. Comparison with previous studies

Here we compare our results with the studies of BC10 and HOK11. We use two particular examples: the long-dashed line in Fig. 2 of BC10 (episodic and cold accretion) and the case of “mC5-C” in HOK11 (steady and cold accretion). The other settings are summarized in Table 1. Quantitatively, their results differ: HOK11 obtain radii that are often about two times larger than those of BC10. For example, in HOK11, the radius when M = 0.9 M is ~ 1.3 R (see their Fig. 2) while for the same mass, it is ~ 0.6 R in BC10.

An important difference between these calculations is that the assumed deuterium contents differ by a factor 1.75. As shown in the previous section, this difference has a large impact on the evolution (see Fig. 3). Figure 6 compares the radius evolution in BC10 and HOK11 to our calculations with the same hypotheses. Given the fact that different stellar models are used (including different values of αMLT – see Table 1), there is excellent agreement between our calculations and those of BC10 and HOK11. The difference between the former and the latter is indeed caused by the different assumed deuterium abundances, XD = 20 ppm and 35 ppm, respectively. The figure also compares two calculations with episodic and steady accretion for XD = 20 ppm. The difference between these two calculations is much smaller than between calculations with different deuterium abundances: the evolution of the PMS star is principally governed by the mean accretion rate and by the deuterium abundance of accreted material and much less by the nature (steady or episodic) of the accretion.

4. Hot and warm accretion

In this section we explore how heat injection by accretion affects stellar evolution. First, we assume that accretion energy is redistributed uniformly in the star and examine variations with the accretion efficiency parameter ξ defined in Eq. (1). Then, we examine, with the parameter mke (see Eq. (3)), how assumptions concerning where the accretion energy is released affect the PMS evolution tracks. In this section the parameters other than ξ and mke are set to their default values (see Table 1).

4.1. Early evolution with a varying ξ

Figure 7 shows the evolution tracks obtained with the heat injection parameter in the range ξ = 0–0.5. When comparing simulations with ξ = 0.5 and ξ = 0, the radii differ by 1 order of magnitude and the luminosities differ by 2 orders of magnitude. These differences in evolutionary tracks are larger and concern a wider range of effective temperatures than those due to a different deuterium content. Thus, the most important parameter controlling the PMS evolution is the amount of heat injected by accretion.

thumbnail Fig. 6

Cold-accretion (ξ = 0) evolution of the stellar radius with time for the studies of HOK11 (dashed green) and BC10 (dashed red) and for our results with steady accretion and XD = 35 ppm (plain green), XD = 20 ppm (dotted red). A model with episodic accretion (as described in Sect. 3.2) and XD = 20 ppm is also shown (plain red).

The evolutionary track of accreting stars corresponds to the “birthline” proposed by Stahler (1983). As BVC12 pointed out, Fig. 7b shows that the birthlines strongly depend on heat injection ξ, implying that the concept of a definitive birthline may be elusive. The radius when the accretion is completed (i.e., at 105 yr) in the case of ξ = 0.5 is almost the same as obtained by Stahler & Palla (2005) for spherical accretion.

thumbnail Fig. 7

Top panel: radius evolution for various heat injection efficiencies, ξ = 0.5 (solid line), 0.1 (dashed), 0.05 (dotted), 0.01 (dot-dashed), and 0 (double dot-dashed). The classical PMS evolution, which is indicated by the thin solid line and starts at the asterisks, has a substantial overlap with the lines of ξ = 0.5. Bottom panel: evolutionary tracks with different ξ ranging 0.5 (red line), 0.1 (green), 0.05 (blue), 0.01 (magenta) and 0 (cyan). The solid, dotted, and dashed lines represent the cases with uniform heat injection, mke = 1, and mke = 0.1, respectively. In both panels, we adopt XD = 20 ppm and Rini = 1.5 R.

Figure 7a shows that before deuterium burning (at 104 to 105 yr), the radius initially increases for ξ = 0.5, decreases slightly for ξ = 0.1, and shows a pronounced decrease for ξ ≲ 0.05. Again, this behavior can be understood using the energy equation as in Appendix B. The situation corresponds to a case when nuclear reactions and radiative cooling may be neglected (i.e., Lnuc,LLacc). The energy equation may thus be shown to yield R=(21ξC)M,\begin{equation} \frac{\dot{R}}{R_\star} = \left( 2-\frac{1-\xi}{C}\right) \frac{\dot{M}}{M_\star}, \label{eq:hot_energy} \end{equation}(6)where C is a constant as a function of the specific heat ratio and the polytropic index (see Appendix C). Therefore, the radius should be constant if ξ = 1−2C. In the case of a fully convective star with monoatomic gas (C = 3/7), the critical ξ for the constant radius is 1/7.

In a second phase, for ξ ≲ 0.1, deuterium burning sets in and leads to a limited increase of the stellar radius. However, this does not occur in the case of ξ = 0.5. First, in this case, deuterium fusion starts only at ~ 1.2 × 105 yr, i.e., after the accretion is completed because of the low central temperature during the accreting phase due to the large radius. Moreover, since Lnuc does not exceed L, deuterium burning only delays the contraction from radiative cooling, and the star does not expand (see Appendix C and Eq. (C.7)).

4.2. Importance of the location of heat injection

So far, we assumed a uniform distribution of injected heat as in BC10 (see Eq. (2)). However, this assumption may not be valid, especially in stars with a large radiative core. Moreover, as described in Sect. 2.6, recent radiation-hydrodynamic simulations by Geroux et al. (2016) showed that the uniform heat redistribution may not be accurate. We now examine the case in which the accretion heat is injected only in surface regions (see Eq. (3)).

Figure 7b shows that the evolution of an accreting solar-mass star strongly depends on the assumed location of the heat injection region. In the case of ξ = 0.05 (blue lines), the luminosity with mke = 0.1 is up to about one order of magnitude larger than the case of the uniform distribution and is almost the same as the case with ξ = 0.3 and uniform distribution. This dependence mainly comes from the assumption that the entropy of accreting materials are the same as the stellar surface (see Sect. 2.6). If mke is small and ξ> 0, the accretion heat is injected only in the outer envelope. This causes εadd to be large at these locations, which, according to Eq. (4), leads to an increase of the stellar luminosity and hence of the specific entropy in this region. The surface entropy then becomes larger than in the case of a large mke, or equivalently, of a uniform injection of accretion heat. The assumption that any added mass has the same entropy as the stellar photosphere effectively leads to accreting material with a higher entropy and therefore to a larger radius.

The assumptions of a uniform or linear injection of accretion energy and of continuity of entropy between the photosphere and accreted material are questionable. However, we can see that the uniform model with low ξ effectively represents one low extreme model. For high ξ values (ξ ≳ 0.3), we expect evolution models to evolve very rapidly initially, losing memory of the initial conditions and resembling standard evolutionary tracks. In that sense, although radiation hydrodynamic simulations would be desirable, the uniform model defined by Eq. (2)is a useful simplification to represent possible evolutionary tracks.

5. Implication for the evolutionary tracks in the H-R diagram

In this section we compare our results with observations in the H-R diagrams. We now turn to isochrones by integrating the evolutions of various final masses as shown in Sect. 3.5. We give special attention to whether these new evolutionary tracks can explain the luminosity spread problem for young stellar objects (YSOs), i.e., the fact that for a given cluster, stars are spread over a relatively wide range of luminosities instead of forming a well-defined luminosity-effective temperature relation as would be expected for stars of similar ages and compositions. Although the consequences of accretion in the H-R diagrams have been investigated by previous works (BCG09, HOK11 and BVC12), we choose here to attempt to constrain the values of ξ that are in agreement with the observations of young clusters.

5.1. The PMS luminosity spread problem

5.1.1. Observational constraints

The luminosity spread of PMS stars has been a matter of debate for decades (see the review of Hillenbrand 2009; Jeffries 2012; Soderblom et al. 2014). This spread is seen almost ubiquitously in star-forming regions and young clusters. Three types of explanations have been proposed: (i) observational or astrophysical uncertainties (e.g., Hartmann 2001); (ii) an intrinsic age spread (e.g., Palla & Stahler 2000; Inutsuka et al. 2015) and (iii) physical processes (e.g., Chabrier et al. 2007; Baraffe et al. 2009). Determining the reason for this spread is important for our understanding of star formation.

It is difficult to determine the luminosity of young PMS stars because it is subject to the observational (e.g., differential extinction, reddening, distance, and cluster membership) and astrophysical (e.g., circumstellar material and its accretion, unresolved binary and variability) uncertainties (e.g., Hartmann 2001). However, the contribution to the luminosity spread by each uncertainty has been quantitatively estimated (e.g., Reggiani et al. 2011; Burningham et al. 2005) and various authors claimed that the sum of their contributions is smaller than the observed luminosity spread. Moreover, Jeffries (2007) found that the projected radii, which are less affected by observational uncertainties, instead of luminosities, also spread widely. These results suggest that the luminosity spread is genuine.

A luminosity spread of 0.2–0.3 dex would correspond to a ~ 0.4 dex spread in the ages of PMS stars. This could be explained theoretically; for example, Inutsuka et al. (2015) propose that stars are formed by the recurrent compressions by expanding bubbles and that consequently the members of a cluster are not necessarily formed in a short period of time.

Instead, cold accretion leads to much smaller stellar radii and luminosities than the classical, non-accreting models and this can explain the luminosity spread for stars Teff ≳ 3500 K (see BCG09; HOK11; BVC12; see also Sect. 3.1). For stars of lower temperatures, the results depend on the size of the assumed seed radius and differ between the different studies (see Sect. 3.4).

After showing how the PMS isochrones depend on ξ, XD, and Rini, we try to estimate possible values of ξ that are compatible with the observations. In this section, we assume a uniform injection of accretion heat (see Sect. 5.1.5).

5.1.2. Isochrones as a function of assumed ξ

In Fig. 8 we compare the PMS stars observed in several young clusters to our theoretical isochrones for the cases ξ = 0, 0.05, 0.1, and 0.5. The remaining parameters are our fiducial values (XD = 20 ppm, Rini = 1.5 R, and Mini = 0.01 M). We use the following observational data of young stars: ρ Ophiuchus (Gatti et al. 2006), σ Orionis (Gatti et al. 2008), Taurus and Chamaeleon I (Muzerolle et al. 2005), Taurus-Auriga (Kenyon & Hartmann 1995), and Orion nebula cluster (ONC; Da Rio et al. 2010).

thumbnail Fig. 8

Comparison of 0.3 (solid line), 1 (dashed), 3 (dotted), and 10 (dot-dashed) million-year isochrones for the cases ξ = 0.5, 0.1, 0.05, and 0 from top to bottom in the H-R diagram with observed PMS stars: ONC (red, small diamonds), ρ Oph (red squares), σ Ori (green circles), Tau-Aur (blue triangles), and Tau and Cha I (blue, inverted triangles). Each isochrone is comprised of evolutionary tracks of 0.01, 0.25, 0.05, 0.1, 0.3, 0.5, 0.8, and 1 M stars denoted by pluses, crosses, asterisks, open squares, open circles, open triangles, open inverted triangles, and open tilted squares, respectively. The initial radius is 1.5 R and XD = 20 ppm. The isochrones with varying ξ cover the wide range of the region where the observed stars exist in the high-temperature regime.

The first panel shows that classical evolutionary tracks with ξ = 0.5 would require a large spread of ages to explain all clusters. The ages of stars would need to range between 0.3 and 10 Myr in ρ Oph, ONC, and Tau-Aur and between 1 and 10 Myr in the other clusters to reproduce the observed luminosity spread. Even in that case, a few stars are underluminous and can be explained only by invoking that their luminosity has been underestimated. Gatti et al. (2006) estimated that the errors in the luminosities of these stars is ~ 0.5 dex, which indeed means that this is a possibility.

Other fixed values of ξ do not allow us to find a better solution. When ξ = 0.1, the most luminous stars cannot be explained anymore. When ξ ≤ 0.05 the slope of the isochrones becomes inconsistent with the ensemble of observational data points. At the same time, the ξ = 0 isochrones are characterized by very low luminosities and can explain the stars with the lowest luminosities observed in ρ Oph.

Conversely, if one assumes that stars within a cluster are coeval, a distribution of the values of ξ within a cluster can be invoked to explain the luminosity spread, as proposed by BCG09, HOK11, and BVC12. However, this approach fails to reproduce the luminosity spread of very-low-mass stars in the Tau and Cha I. We revisit this in Sect. 5.1.4.

5.1.3. Dependence on deuterium content

As described in Sect. 3.3, even the deuterium fraction of the present-day local ISM is still under debate. We now investigate how isochrones depend on deuterium content.

Figure 9 shows 1 Myr isochrones obtained with different values of ξ and with either XD = 20 ppm or 35 ppm. The isochrones obtained ξ = 0.5 are almost independent of XD. However, isochrones obtained for ξ = 0 differ very significantly when XD changes. As described in Sect. 3, this is because in the case of cold accretion, deuterium burning regulates the PMS evolution. But Fig. 9 also shows that the isochrones remain almost parallel. For example, the isochrone with (ξ,XD) = (0.1,20 ppm) is very similar to that with (0,35 ppm). This is because both ξ and XD control the specific entropy of the accreted material.

In the context of a variable ξ value, a low abundance of deuterium yields a more important spread of solutions than a high abundance. We do not expect XD to vary within a cloud because of turbulent mixing. However, if XD varies from one cluster to the next, this would yield a more important luminosity spread for clusters with low XD values.

thumbnail Fig. 9

Comparison of 1 Myr isochrones for the cases XD = 35 (red lines) and 20 ppm (blue) ranging ξ from 0 (dotted) to 0.5 (solid), respectively. The initial radius is 1.5 R. The points for observation and isochrones are the same as in Fig. 8. The isochrone with (XD) = (35 ppm,0) is almost the same as that with (20 ppm, 0.1), implying that the different deuterium content has a large impact on the PMS evolution. The thin solid line is a 1 Myr non-accreting isochrone (Baraffe et al. 1998).

5.1.4. Dependence on initial conditions

As described in Sect. 3.4 in the context of cold accretion, the PMS evolutionary tracks depend on the initial conditions, and in particular on the physical characteristics of the initial stellar seed. In Fig. 10, we show the isochrones obtained with Rini of 0.25 R and 3 R, for different values of ξ. The comparison of the results with the same ξ and different Rini shows that the dependence on Rini is significant across the entire effective temperature range, except for hot accretion (ξ ~ 0.5). The dependence is most important for cold accretion (see dotted lines in Fig. 10) and corresponds to the situation described in Sect. 3.4. For intermediate values ξ ~ 0.1 (dotted lines), the dependence on Rini is smaller but still significant. It is most pronounced at small ages (≲ 1 Myr, which corresponds to the K-H timescale) and for low-mass stars; since the total energy of high-mass stars is dominated by the injected energy, Eadd, the initial energy does not matter.

thumbnail Fig. 10

Comparison of isochrones at 0.3, 1, 3, and 10 Myr from top to bottom. Rini is either 3 R (gray lines) or 0.25 R (black). The solid, dashed, and dotted lines represent isochrones with ξ = 0.5,0.1, and 0, respectively. Here XD = 20 ppm. The points are the same as in Fig. 8.

A high value of Rini ≳ 1.5 R is required to explain the most luminous stars, particularly for the low Teff < 3000 K stars. However, invoking a variable Rini including low values Rini ~ 0.25 R is required to explain the observed luminosity spread in the Tau and Cha I clusters by assuming coeval stars.

5.1.5. Resulting constraints

The observed luminosities and effective temperatures of young clusters can hence be explained in the context of coeval stars by assuming that ξ and Rini can vary within a given cluster. This interpretation suggests that ρ Oph, ONC, and Tau-Aur are extremely young (~ 0.3 Myr); σ Ori is slightly older at ~ 1 Myr5. However, inaccurate estimations of stellar luminosities, membership issues, and other problems can affect these age estimates.

We can see in Sect. 5.1.2 and Fig. 8 that for a majority of stars ξ ≳ 0.1 for XD = 20 ppm and Rini = 1.5 R considering the following two facts: the slopes of the isochrones are not compatible with the observations for ξ ≤ 0.05 and the number of underluminous stars is small in Teff = 3500–4500 K. For higher values of XD, this constraint on ξ should be relaxed somewhat. For lower values of Rini, we would have to impose ξ to be even closer to 0.5. Even with higher values of Rini, ξ ≲ 0.1 would still be rare because the slope of the isochrone does not match the observation. It is thus premature at this point and without a proper modeling of where the accretion heat is deposited to attempt constraining ξ more precisely. However, one important conclusion that we can draw is that the evolution tracks cannot deviate too significantly from the standard model: stars with entropies that are too low are ruled out by the observation of young clusters. In that sense, the limit set by the model characterized by a uniform accretion heat redistribution, ξ = 0.1, XD = 20 ppm and Rini = 1.5 R, is useful when considering the range of possibilities in agreement with observational data.

Finally, if indeed star growth is characterized by an inefficient burial of accretion heat (ξ< 0.5), then we should expect that clusters characterized by a lower deuterium mixing ratio XD should have a larger luminosity spread than clusters with a higher XD. This may thus become testable, depending on the possibility to reliably determine D/H ratios in clusters (see Sect. 3.3).

5.2. The case of CoRoT 223992193

We now consider the case of a specific system, the eclipsing binary CoRoT 223992193 (Gillen et al. 2014; Stassun et al. 2014). We expect both stars to have the same age. Therefore, the system provides important constraints to test our evolution models and retrieve ξ independent of the cluster results discussed previously.

Figure 11 shows the constraints on luminosity and effective temperature for both components of the system and compares them to evolutionary tracks for classical (non-accreting) models, and for our models with values of ξ between 0 and 0.1, XD of 20 and 35 ppm, and Rini = 1.5 R. The two stars have masses between 0.6 M and 0.8 M and an age less than 5 Myr. Although the evolutionary tracks are strongly sensitive to the assumed values of ξ, the isochrones remain mostly parallel to each other, except for extremely low values of ξ with small XD. This implies that varying ξ essentially results in a shift in age, where stars with a given effective temperature and luminosity are younger for lower values of ξ.

thumbnail Fig. 11

Evolutionary tracks whose final masses are 0.5, 0.8, and 1 M for the cases of ξ = 0 (green lines), 0.05 (blue), and 0.1 (magenta dot-dashed). The initial radius is 1.5 R and XD = 20 ppm (solid lines) and 35 ppm (dashed). Some are merged to the Hayashi track of the non-accreting evolutionary tracks (black dotted), and others are not. Thus, the mass determination of stars on the Hayashi phase is not affected by the variation in accretion heat. As already discussed, the age determination is affected by ξ. Pluses, crosses, and asterisks denote 1, 5, and 10 million years, respectively. The filled and open circles represent the primary and secondary of the eclipsing binaries CoRoT 223992193, respectively (Stassun et al. 2014). For the cases of XD = 20 ppm, their age is estimated to 1 or 5 million years if ξ = 0.05 or 0.1, respectively.

Quantitatively, the constraints that we can derive on ξ are very similar to those obtained for young clusters (and with the same caveats). With the assumptions of Rini and the uniform distribution, for XD = 20 ppm, we obtain that ξ ≥ 0.05, with the limiting case corresponding to both stars being close to their birthline. For XD = 35 ppm, the initial entropy is larger so that no constraint on ξ can be obtained.

6. Conclusions

The PMS evolution of stars has long been considered a relatively simple theoretical problem governed by the quasi-static contraction of an almost isentropic star. We have seen that this evolution phase may in fact be strongly altered by the fact that the stellar envelope must be accreted onto a protostellar seed and that a significant fraction of the energy may be lost in the accretion shock connected to this process. The goal of this paper was to understand what controls the PMS and to derive constraints from the observations. In order to do so, we used simulations with the stellar evolution code MESA that account for a progressive accretion of material onto a forming star.

We have first shown that beyond classical parameters, such as mass and metallicity, the evolution on the PMS is controlled essentially by three parameters: ξ, the efficiency at which the gravitational energy of the accreted material is transformed into internal energy of the star, XD, the mass mixing ratio of deuterium, and the entropy of the initial stellar seed (or equivalently in the present work its radius Rini for a mass of 0.01 M). The parameter ξ = 0.5 corresponds to the classical models and lead to the formation of stars with a large radius and entropy and an evolution that is essentially independent of the two other parameters, XD and Rini. Progressively lower values of ξ yield stars with (much) smaller radii and a richer ensemble of possibilities in terms of their evolution. In particular, because the entropy of the accreted material is then effectively smaller, the abundance of deuterium in that material becomes very important in deciding whether stars will expand significantly (high XD values) or will retain a small size all the way to the main sequence (for XD ≲ 30 ppm). We showed that the differences between the results of Baraffe et al. (2009) and Hosokawa et al. (2011) can be accurately reproduced and the results from different choices for XD.

We compared the evolutionary tracks to the observations of several young clusters. We confirmed that the spread in luminosities in each cluster could be explained without invoking an age spread but by instead assuming that ξ could vary from one star to the next. A variation of Rini (i.e., stellar entropy at 0.01 M) is also needed to explain underluminous, cool stars in Tau and Cha I. However, the observations indicate that most stars cannot be too low in entropy when they form. This implies that within the uniform accretion model, we can rule out as unlikely those scenarios with low ξ and XD values. Specifically, the model with XD = 20 ppm, Rini = 1.5 R, and ξ = 0.1 sets a useful boundary: stars can have lower entropies only in relatively rare cases. This means that for a majority of stars, stellar evolution cannot differ from the classical evolution tracks beyond the limit set by this limiting model. On the other hand, because of the multi-parameter dependence, we cannot derive an independent constraint on ξ. For example, a model with XD = 35 ppm and ξ = 0 is equivalent in terms of entropy and luminosities in the H-R diagram to the XD = 20 ppm and ξ = 0.1. We found these constraints to be compatible with the observational constraints from the PMS eclipsing binary CoRoT 223992193. One caveat is that if a significant number of stars in the clusters are affected by observational errors (such as L, membership), our constraints is changed.

Separately, we found that the possibility of reliably measuring deuterium abundances in clusters would allow testing for inefficient accretion (i.e., ξ< 0.5). If this is the case, we would expect the spread in luminosity to be larger in clusters with a lower deuterium to hydrogen ratio. Present observations indicate that σ Ori may have a smaller luminosity spread than other clusters. Although many other factors are to be considered, it is possible that this is due to a higher D/H ratio in that cluster.

The PMS evolution of stars is not as simple as once thought and merits further investigation. On the observational side, further insight would be gained through the discovery of more very young eclipsing binaries, the determination of deuterium abundances in clusters, and further observations of young accreting stars. On the theoretical side, the main uncertainties in our calculations are due to extremely simplified outer boundary conditions and the ad hoc prescription used to relate the gravitational energy of the accreted material to the internal energy in the star. Three-dimensional radiation hydrodynamic simulations of a collapsing molecular cloud core with sufficient resolution to resolve the central stellar seed are needed.

A consequence of this more complex PMS evolution is that the stellar interior is not necessarily fully convective during most of this phase. This may have strong implications to understand the chemical composition of stars and connect measurements of stellar compositions to the formation of planets. We will investigate this issue in our next paper.


1

Our ξ is equivalent to αϵ in Baraffe et al. (2009, 2012).

2

These values are estimated using the number ratio of the deuterium over hydrogen, (D/H). Since the hydrogen mass fraction X ≃ 0.7 and the mass ratio of deuterium to hydrogen is about two, XD ≃ 1.4(D/H).

3

Two methods may be used to constrain the protosolar deuterium abundance: using either the deuterium content in the Jovian atmosphere (Lellouch et al. 2001; Guillot 1999) or the enhanced 3He measured in the solar wind and compared to the Jovian atmosphere (Heber et al. 2008).

4

The luminosity in the present paper is always the stellar intrinsic luminosity and does not include the accretion luminosity emitted from the shock surface.

5

We do not derive the age of Tau and Cha I because the data in Muzerolle et al. (2005) is an assembly of young stars in two star-forming regions.

6

Using three-dimensional atmosphere model, the metallicity in the solar atmosphere is estimated to be ~ 0.0134, which is reduced from the classical value (~ 0.02) by about 70% (see Asplund et al. 2009).

7

The total luminosity of protostars is the sum of the intrinsic luminosity, L, and the radiation from the accretion shock front (i.e., LaccLadd).

8

In addition, another derivation is possible using a polytropic analysis (Chandrasekhar 1967).

9

This is equivalent to the condition LLacc as shown in Fig. B.1.

10

Vaytet et al. (2013) indicate that the difference probably results from different opacity tables.

Acknowledgments

We express our gratitude to S. Inutsuka, T. Hosokawa, P. Morel, T. Nakamoto, M. Kuzuhara, and M. Ikoma for fruitful discussions and comments. Bill Paxton and Dean Townsley kindly helped M.K. use the stellar-evolution code MESA. We appreciate the critical and constructive comments of the referees, which helped us to improve this paper. M.K. is supported by Grant-in-Aid for JSPS Fellows Grant Number 24·9296, MEXT of Japan (Grant: 23244027) and Foundation for Promotion of Astronomy.

References

  1. Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., & Schweitzer, A. 2001, ApJ, 556, 357 [NASA ADS] [CrossRef] [Google Scholar]
  2. Amelin, Y., Krot, A. N., Hutcheon, I. D., & Ulyanov, A. A. 2002, Science, 297, 1678 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  3. Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nucl. Phys. A, 656, 3 [NASA ADS] [CrossRef] [Google Scholar]
  4. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  5. Bahcall, J. N., Basu, S., Pinsonneault, M., & Serenelli, A. M. 2005, ApJ, 618, 1049 [NASA ADS] [CrossRef] [Google Scholar]
  6. Baraffe, I., & Chabrier, G. 2010, A&A, 521, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  7. Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 337, 403 [NASA ADS] [Google Scholar]
  8. Baraffe, I., Chabrier, G., & Gallardo, J. 2009, ApJ, 702, L27 [NASA ADS] [CrossRef] [Google Scholar]
  9. Baraffe, I., Vorobyov, E., & Chabrier, G. 2012, ApJ, 756, 118 [NASA ADS] [CrossRef] [Google Scholar]
  10. Basu, S., & Antia, H. M. 2004, ApJ, 606, L85 [NASA ADS] [CrossRef] [Google Scholar]
  11. Burningham, B., Naylor, T., Littlefair, S. P., & Jeffries, R. D. 2005, MNRAS, 363, 1389 [NASA ADS] [CrossRef] [Google Scholar]
  12. Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., & Salaris, M. 2007, ApJ, 661, 1094 [NASA ADS] [CrossRef] [Google Scholar]
  13. Castelli, F., & Kurucz, R. L. 2004, IAU Symp. 210, poster A20 [Google Scholar]
  14. Chabrier, G., Gallardo, J., & Baraffe, I. 2007, A&A, 472, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Chambers, J. E. 2010, ApJ, 724, 92 [NASA ADS] [CrossRef] [Google Scholar]
  16. Chandrasekhar, S. 1967, An introduction to the study of stellar structure (New York: Dover Publications) [Google Scholar]
  17. Cox, J., & Giuli, R. 1968, Principles of Stellar Structure (New York: Gordon and Breach), 401 [Google Scholar]
  18. Da Rio, N., Robberto, M., Soderblom, D. R., et al. 2010, ApJ, 722, 1092 [NASA ADS] [CrossRef] [Google Scholar]
  19. Dunham, M. M., & Vorobyov, E. I. 2012, ApJ, 747, 52 [NASA ADS] [CrossRef] [Google Scholar]
  20. Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, ApJ, 623, 585 [NASA ADS] [CrossRef] [Google Scholar]
  21. Gatti, T., Testi, L., Natta, A., Randich, S., & Muzerolle, J. 2006, A&A, 460, 547 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Gatti, T., Natta, A., Randich, S., Testi, L., & Sacco, G. 2008, A&A, 481, 423 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  23. Geroux, C., Baraffe, I., Viallet, M., et al. 2016, A&A, 588, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  24. Gillen, E., Aigrain, S., McQuillan, A., et al. 2014, A&A, 562, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  25. Grevesse, N., & Sauval, A. J. 1998, Space Sci. Rev., 85, 161 [NASA ADS] [CrossRef] [Google Scholar]
  26. Guillot, T. 1999, Planet. Space Sci., 47, 1183 [NASA ADS] [CrossRef] [Google Scholar]
  27. Guillot, T., Ida, S., & Ormel, C. W. 2014, A&A, 572, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  28. Hartmann, L. 2001, AJ, 121, 1030 [NASA ADS] [CrossRef] [Google Scholar]
  29. Hartmann, L., Cassen, P., & Kenyon, S. J. 1997, ApJ, 475, 770 [NASA ADS] [CrossRef] [Google Scholar]
  30. Hauschildt, P. H., Allard, F., & Baron, E. 1999a, ApJ, 512, 377 [NASA ADS] [CrossRef] [Google Scholar]
  31. Hauschildt, P. H., Allard, F., Ferguson, J., Baron, E., & Alexander, D. R. 1999b, ApJ, 525, 871 [NASA ADS] [CrossRef] [Google Scholar]
  32. Hayashi, C. 1961, PASJ, 13, 450 [NASA ADS] [Google Scholar]
  33. Heber, V. S., Baur, H., Bochsler, P., et al. 2008, in Lunar and Planetary Inst. Tech. Rep., Lunar and Planetary Science Conference, 39, 1779 [NASA ADS] [Google Scholar]
  34. Hébrard, G., Tripp, T. M., Chayer, P., et al. 2005, ApJ, 635, 1136 [NASA ADS] [CrossRef] [Google Scholar]
  35. Henyey, L., Vardya, M. S., & Bodenheimer, P. 1965, ApJ, 142, 841 [NASA ADS] [CrossRef] [Google Scholar]
  36. Herwig, F. 2000, A&A, 360, 952 [NASA ADS] [Google Scholar]
  37. Hillenbrand, L. A. 2009, in IAU Symp. 258, eds. E. E. Mamajek, D. R. Soderblom, & R. F. G. Wyse, 81 [Google Scholar]
  38. Hosokawa, T., & Omukai, K. 2009, ApJ, 691, 823 [NASA ADS] [CrossRef] [Google Scholar]
  39. Hosokawa, T., Offner, S. S. R., & Krumholz, M. R. 2011, ApJ, 738, 140 [NASA ADS] [CrossRef] [Google Scholar]
  40. Hueso, R., & Guillot, T. 2005, A&A, 442, 703 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  41. Inutsuka, S.-i. 2012, Prog. Theor. Exp. Phys., 2012, 010000 [Google Scholar]
  42. Inutsuka, S.-i., Machida, M. N., & Matsumoto, T. 2010, ApJ, 718, L58 [NASA ADS] [CrossRef] [Google Scholar]
  43. Inutsuka, S.-i., Inoue, T., Iwasaki, K., & Hosokawa, T. 2015, A&A, 580, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  44. Jeffries, R. D. 2007, MNRAS, 381, 1169 [NASA ADS] [CrossRef] [Google Scholar]
  45. Jeffries, R. D. 2012, Are There Age Spreads in Star Forming Regions?, eds. A. Moitinho, & J. Alves (Berlin, Heidelberg: Springer-Verlag), 163 [Google Scholar]
  46. Kenyon, S. J., & Hartmann, L. 1995, ApJS, 101, 117 [NASA ADS] [CrossRef] [Google Scholar]
  47. Larson, R. B. 1969, MNRAS, 145, 271 [NASA ADS] [CrossRef] [Google Scholar]
  48. Lellouch, E., Bézard, B., Fouchet, T., et al. 2001, A&A, 370, 610 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  49. Linsky, J. L., Draine, B. T., Moos, H. W., et al. 2006, ApJ, 647, 1106 [NASA ADS] [CrossRef] [Google Scholar]
  50. Mahaffy, P. R., Donahue, T. M., Atreya, S. K., Owen, T. C., & Niemann, H. B. 1998, Space Sci. Rev., 84, 251 [NASA ADS] [CrossRef] [Google Scholar]
  51. Martin, R. G., Lubow, S. H., Livio, M., & Pringle, J. E. 2012, MNRAS, 423, 2718 [NASA ADS] [CrossRef] [Google Scholar]
  52. Masunaga, H., & Inutsuka, S.-i. 2000, ApJ, 531, 350 [NASA ADS] [CrossRef] [Google Scholar]
  53. Meléndez, J., Asplund, M., Gustafsson, B., & Yong, D. 2009, ApJ, 704, L66 [NASA ADS] [CrossRef] [Google Scholar]
  54. Mercer-Smith, J. A., Cameron, A. G. W., & Epstein, R. I. 1984, ApJ, 279, 363 [NASA ADS] [CrossRef] [Google Scholar]
  55. Muzerolle, J., Luhman, K. L., Briceño, C., Hartmann, L., & Calvet, N. 2005, ApJ, 625, 906 [NASA ADS] [CrossRef] [Google Scholar]
  56. Nelder, J. A., & Mead, R. 1965, Comput. J., 7, 308 [CrossRef] [Google Scholar]
  57. Palla, F., & Stahler, S. W. 1992, ApJ, 392, 667 [NASA ADS] [CrossRef] [Google Scholar]
  58. Palla, F., & Stahler, S. W. 2000, ApJ, 540, 255 [NASA ADS] [CrossRef] [Google Scholar]
  59. Paxton, B., Bildsten, L., Dotter, A., et al. 2011, ApJS, 192, 3 [Google Scholar]
  60. Paxton, B., Cantiello, M., Arras, P., et al. 2013, ApJS, 208, 4 [NASA ADS] [CrossRef] [Google Scholar]
  61. Paxton, B., Marchant, P., Schwab, J., et al. 2015, ApJS, 220, 15 [Google Scholar]
  62. Prantzos, N. 2007, Space Sci. Rev., 130, 27 [NASA ADS] [CrossRef] [Google Scholar]
  63. Ramírez, I., Meléndez, J., & Asplund, M. 2009, A&A, 508, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Ramírez, I., Meléndez, J., Cornejo, D., Roederer, I. U., & Fish, J. R. 2011, ApJ, 740, 76 [NASA ADS] [CrossRef] [Google Scholar]
  65. Reggiani, M., Robberto, M., Da Rio, N., et al. 2011, A&A, 534, A83 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Rogers, F. J., & Nayfonov, A. 2002, ApJ, 576, 1064 [Google Scholar]
  67. Saumon, D., Chabrier, G., & van Horn, H. M. 1995, ApJS, 99, 713 [NASA ADS] [CrossRef] [Google Scholar]
  68. Seaton, M. J. 2005, MNRAS, 362, L1 [NASA ADS] [CrossRef] [Google Scholar]
  69. Shu, F. H. 1977, ApJ, 214, 488 [NASA ADS] [CrossRef] [Google Scholar]
  70. Soderblom, D. R., Hillenbrand, L. A., Jeffries, R. D., Mamajek, E. E., & Naylor, T. 2014, in Protostars and Planets VI (Tucson: University of Arizona Press), 219 [Google Scholar]
  71. Stahler, S. W. 1983, ApJ, 274, 822 [NASA ADS] [CrossRef] [Google Scholar]
  72. Stahler, S. W. 1988, ApJ, 332, 804 [NASA ADS] [CrossRef] [Google Scholar]
  73. Stahler, S. W., & Palla, F. 2005, The Formation of Stars (Wiley-VCH) [Google Scholar]
  74. Stahler, S. W., Shu, F. H., & Taam, R. E. 1980, ApJ, 241, 637 [NASA ADS] [CrossRef] [Google Scholar]
  75. Stahler, S. W., Palla, F., & Salpeter, E. E. 1986, ApJ, 302, 590 [NASA ADS] [CrossRef] [Google Scholar]
  76. Stamatellos, D., Whitworth, A. P., Bisbas, T., & Goodwin, S. 2007, A&A, 475, 37 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Stassun, K. G., Feiden, G. A., & Torres, G. 2014, New Astron. Rev., 60, 1 [NASA ADS] [CrossRef] [Google Scholar]
  78. Steigman, G. 2006, Int. J. Mod. Phys. E, 15, 1 [Google Scholar]
  79. Sugimoto, D., & Nomoto, K. 1975, PASJ, 27, 197 [NASA ADS] [Google Scholar]
  80. Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2013, ApJ, 763, 6 [NASA ADS] [CrossRef] [Google Scholar]
  81. Townsley, D. M., & Bildsten, L. 2004, ApJ, 600, 390 [NASA ADS] [CrossRef] [Google Scholar]
  82. Vaytet, N., Chabrier, G., Audit, E., et al. 2013, A&A, 557, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Vidal-Madjar, A., & Gry, C. 1984, A&A, 138, 285 [NASA ADS] [Google Scholar]
  84. Vorobyov, E. I., & Basu, S. 2005, ApJ, 633, L137 [NASA ADS] [CrossRef] [Google Scholar]
  85. Vorobyov, E. I., & Basu, S. 2010, ApJ, 719, 1896 [NASA ADS] [CrossRef] [Google Scholar]
  86. Winkler, K.-H. A., & Newman, M. J. 1980, ApJ, 236, 201 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: χ2 test for the input parameters

In this study we chose the input parameters that reproduce the observed solar quantities. In addition to the radius and luminosity, the internal structure and surface composition are constrained by helioseismic and spectroscopic analyses (Bahcall et al. 2005; Basu & Antia 2004; Grevesse & Sauval 1998, and references therein). Although the solar metallicity remains a matter of debate6 it is beyond the scope of this paper.

We performed a χ2 test to find the best initial settings using the “Nelder-Mead simplex algorithm” (Nelder & Mead 1965). The input parameters are the initial composition (Xini, Yini, and Zini; see Sect. 2.3), the mixing-length parameter, αMLT, and the overshoot mixing parameter, fov (Herwig 2000; Paxton et al. 2011, see Sect. 2.1).

We calculated the evolution of 1 M stars from PMS phase with changing these parameters to minimize the χ2 value at the solar age, which is assumed to be 4.567 Gyr (Amelin et al. 2002). The results of the χ2 test listed in Table A.1 are used in this paper.

Table A.1

Results of the χ2 test to reproduce solar values.

Appendix B: Evolution and underlying physics in the cold accretion case

In this Appendix, we explain the basic behavior and its underlying physics in the case of cold accretion shown in Sect. 3. The evolution can be split into five phases as shown in Fig. 1. We explain each phase below.

I. The adiabatic contraction phase:

in this phase, the star shrinks while increasing its mass. The radius evolution can be fitted by RM1/3.\appendix \setcounter{section}{2} \begin{equation} R_\star \propto M_\star^{-1/3}. \label{eq:R-M} \end{equation}(B.1)Indeed, we show in Appendix C that this relation is verified for a perfect gas star that accretes gas with the same entropy. In the current settings, the accretion is intense enough that the Kelvin-Helmholtz (K-H) timescale, τKH ≡ | Etot | /L, is much longer than the accretion timescale, τaccM/, where Etot is the total energy of the star and L is the stellar intrinsic luminosity7. This is equivalent to LaccL as shown in Fig. B.1. The accretion luminosity is defined as LaccGM/R,\appendix \setcounter{section}{2} \begin{eqnarray} L\sub{acc} \equiv GM_\star\dot{M}/R_\star, \label{eq:Lacc} \end{eqnarray}(B.2)which is for example 31.3 L in the case of M = 0.1 M, R = 1 R, and = 10-5M/ yr.

thumbnail Fig. B.1

Evolution of the intrinsic luminosity (L, solid line), the energy production rate by thermonuclear reaction (Lnuc, dashed) and the seventh part of the gravitational energy of the accreting material (17Lacc\hbox{$\frac{1}{7}L\sub{acc}$}, dotted). The shaded regions and labels from (I) to (V) are the same as in Fig. 1. During phase (II), Lnuc exceeds 17Lacc\hbox{$\frac{1}{7}L\sub{acc}$} and the star expands. The flat Lnuc during phase (III) results from the steady and instantaneous fusion of the newly accreted deuterium.

II. The deuterium-burning phase:

after the central temperature exceeds ~ 106 K, deuterium fusion affects the evolution (see Sect. 2.2). In the current settings, it happens at t = 7 × 103 yr and M ≃ 0.08 M. As described in Sect. 2.2, the energy production rate of deuterium burning has a strong temperature dependence, that is, εnuc ∝ (T/ 106 K)11.8. This strong temperature sensitivity is responsible for a rapid expansion of the star through the so-called “thermostat effect” (e.g., Stahler 1988).

Indeed, after the ignition of deuterium fusion, the central temperature is maintained constant at about 106 K by the fact that any temperature increase would result in an expansion of the deuterium burning region, and its adiabatic cooling and a temperature decrease would be balanced by adiabatic heating, respectively. With the approximate relations for a perfect gas star in Appendix C, TcPc/ρcM/R and therefore a constant central temperature implies that a mass increase results in an expansion of the star.

However, the central temperature is not exactly constant so we must test that the energy produced by deuterium burning is sufficient to cause the expansion. In Appendix C, we derive Eq. (C.7), which shows how the rate of change of the stellar radius of a fully convective perfect gas star depends on the mass accretion rate, intrinsic luminosity, and nuclear energy production rate. Our derivation follows Hartmann et al. (1997). By neglecting the intrinsic luminosity (since the Kelvin-Helmholtz timescale is long compared to relevant timescales in this phase), R=13M+73RLnucGM2·\appendix \setcounter{section}{2} \begin{equation} \frac{\dot{R}}{R_\star} = -\frac{1}{3}\frac{\dot{M}}{M_\star}+\frac{7}{3}\frac{R_\star L\sub{nuc}}{GM_\star^2}\cdot \label{eq:lnuccrit} \end{equation}(B.3)Therefore, the condition for the expansion is given by Lnuc>17Lacc.\appendix \setcounter{section}{2} \begin{equation} L\sub{nuc} > \frac{1}{7}L\sub{acc}. \label{eq:lnuccrit2} \end{equation}(B.4)In the current phase, since this condition is satisfied owing to the vigorous burning of the pre-existent deuterium, > 0.

III. A second contraction phase:

after ~ 1.5 × 104 yr, pre-existing deuterium has been burned up and deuterium fusion relies on the accretion of fresh deuterium. In this phase, the star contracts again. This is because the mass accretion rate of freshly accreted deuterium, XD, is not sufficient to compensate for the compression by accretion. The maximum energy production rate of deuterium burning, LD, max, is calculated by assuming that it is instantaneously burned, LD,max=8.65L(10-5M/yr)(XD2.0×10-5)·\appendix \setcounter{section}{2} \begin{equation} L\sub{D,max}=8.65~L_\odot \brafracket{\dot{M}}{10^{-5}~\Msun/{\textrm{yr}}}\brafracket{X\sub{D}}{2.0\times10^{-5}}\cdot \end{equation}(B.5)This may be compared to the accretion luminosity Lacc defined in Eq. (B.2), but for a mass that is now M ~ 0.5 M. We thus obtain LD,max/Lacc ~ 0.05 for standard values of the parameters. According to Eq. (B.4), this shows that the burning of accreted deuterium cannot prevent the adiabatic contraction, independent of the accretion rate. In this phase, thermonuclear energy production is dominated by deuterium burning, implying that LnucLD,max until the accretion ceases at 105 yr (see Fig. B.1).

IV. The swelling phase:

after the accretion ceases, the radius remains nearly constant for several million years. This timescale is then determined by what happens in the deep interior and is shorter than the photospheric K-H timescale of ~ 107108 yr. This is because the luminosity in the deep interior is much larger than that at the surface because of the absorption of energy in the stellar interior. Thus, the internal thermal timescale, which is ~GMr2/(2rl)\hbox{${\sim}GM_r^2/(2rl)$}, is down to a few million years in the deep interior, where r and l are the radius and luminosity. For most of this phase, the dominant source of the luminosity is the energy release in each mass shell, i.e., T∂S/∂t in Eq. (4).

After several million years, the star expands. A “luminosity wave” (Stahler et al. 1986) gradually propagates from the interior to the atmosphere and accompanies an increase in temperature and decrease in opacity of the outer layers. The redistribution of entropy in the star due to the luminosity wave changes the internal structure (e.g., the polytropic index n) and then causes the stellar expansion (Hosokawa & Omukai 2009). These profiles are shown in Fig. B.2.

thumbnail Fig. B.2

Profiles of luminosity (top panel), entropy (middle), and radius (bottom) in the stellar interior at 0.1 (red solid line; just after the accretion is completed), 1.0 (green dashed), 2.9 (blue dotted), and 5.1 (magenta dot-dashed) million years. As the luminosity propagates toward the surface, the entropy is redistributed. Crosses represent the base of the surface convective zone; kB and mamu are the Boltzmann constant and the atomic mass unit, respectively.

V. The main sequence:

the star then slightly shrinks from the K-H contraction to enter the main sequence (MS). In this phase, the intrinsic luminosity is almost entirely due to hydrogen burning. The energy equation becomes almost time independent, which means that stars evolve on a much longer timescale. This slow evolution is caused by the change in chemical composition of the central regions due to nuclear energy production. The MS lasts until, in the central regions, hydrogen is exhausted.

Appendix C: Analytical relations

Here we derive the mass-radius relationship of Eq. (B.1) in two ways8. First, we use the characteristic density and pressure and the entropy following Stahler (1988). The characteristic density and pressure of the star in the hydrostatic equilibrium are given by ˜ρM/R3,˜PM2/R4.\appendix \setcounter{section}{3} \begin{eqnarray} \tilde{\rho} \propto M_\star/R_\star^3, ~\tilde{P}&\propto M_\star^2/R_\star^4 \label{eq:rhoP}. \end{eqnarray}(C.1)The entropy is given by S=cVln(Pργad)+S0\hbox{$S=c_V\ln{\left(\frac{P}{\rho^{\gamma\sub{ad}}}\right) }+S_0$}, where cV is the specific heat at constant volume, γad = cP/cV the specific heat ratio, and S0 is the constant. Substituting Eq. (C.1) into this equation, we obtain R=M2γad3γad4exp(13γad4SS0cV)·\appendix \setcounter{section}{3} \begin{equation} R_\star=M_\star^{-\frac{2-\gamma\sub{ad}}{3\gamma\sub{ad}-4}} \exp \left( \frac{1}{3\gamma\sub{ad}-4} \frac{S-S_0}{c_V} \right)\cdot \end{equation}(C.2)In the case of the monatomic ideal gas (γad = 5/3), we obtain R=M1/3exp[23μ(SS0)],\appendix \setcounter{section}{3} \begin{equation} R_\star= M_\star^{-1/3} \exp\left[ \frac{2}{3}\frac{\mu}{\mathcal{R}} (S-S_0) \right], \label{eq:R-M2} \end{equation}(C.3)where μ is the mean molecular weight and is the gas constant. This equation shows that the radius is determined only by entropy and mass. In isentropic stars, RM1/3\hbox{$R_\star\propto M_\star^{-1/3}$} (Eq. (B.1)).

Secondly, Eq. (B.1) can be derived from the energy equation, which also gives Eq. (B.3) following Hartmann et al. (1997). The total energy of a star is given by the sum of the internal and gravitational energy, Etot = Eint + Eg. The energy conservation is expressed as follows: dEtotdt=L+LnucGMR+Ladd.\appendix \setcounter{section}{3} \begin{equation} \frac{{\mathrm{d}}E\sub{tot}}{\mathrm{d}t} =-L_\star +L\sub{nuc} -\frac{GM_\star\dot{M}}{R_\star}+L\sub{add}. \label{eq:dEdt} \end{equation}(C.4)The right-hand side (RHS) shows the radiative cooling, the energy production by thermonuclear reactions, and gravitational and internal energies of accreting materials. From Eq. (1) here we assume Ladd = ξGM/R. From the virial theorem, Etot=(43γad)Eint=3γad43(γad1)Eg.\appendix \setcounter{section}{3} \begin{equation} E\sub{tot}=(4-3\gamma\sub{ad})E\sub{int}=\frac{3\gamma\sub{ad}-4}{3(\gamma\sub{ad}-1)}E\sub{g}. \label{eq:Etot} \end{equation}(C.5)If we assume the polytropic relation that P(ρ) = Kρ1 + 1 /n, Eg=35nGM2R,\appendix \setcounter{section}{3} \begin{equation} E\sub{g}=-\frac{3}{5-n}\frac{GM_\star^2}{R_\star}, \end{equation}(C.6)where n is the polytropic index. If we define C ≡ (3γad−4)/(γad−1)(5−n), Etot=CGM2/R\hbox{$E\sub{tot}=-CGM_\star^2/R_\star$}. Thus, the total energy evolution in Eq. (C.4) is transformed as R=(21ξC)MRLCGM2+RLnucCGM2·\appendix \setcounter{section}{3} \begin{equation} \frac{\dot{R}}{R_\star} = \left( 2-\frac{1-\xi}{C}\right) \frac{\dot{M}}{M_\star} -\frac{RL_\star}{CGM^2}+\frac{RL\sub{nuc}}{CGM^2}\cdot \label{eq:KHMacc} \end{equation}(C.7)The second term of RHS of Eq. (C.7) corresponds to the inverse of the K-H timescale τKH, which is defined as the typical timescale of the radiative cooling (i.e., τKH ≡ | Etot | /L). During the main-accretion phase, τKH is in general much longer than the accretion timescale τaccM/9. Therefore, the second term of RHS of Eq. (C.7) can be neglected if the first term of RHS is not zero, i.e., 1−ξ ≠ 2C.

In fully convective stars, which consist of the monatomic ideal gas, n = 3/2 and γad = 5/3 and then C = 3/7. In addition, in the case that ξ = 0, we obtain Eq. (B.3). Moreover if LnucLacc, it is transferred as Eq. (B.1).

Appendix D: Dependence on the initial stellar seed mass

In this paper we chose 10-2M as the fiducial value of the initial stellar seed mass. We stress that this value is higher than the second-core mass in recent work of the hydrodynamic collapse of molecular clouds (~ MJup; see Masunaga & Inutsuka 2000; Stamatellos et al. 2007; Tomida et al. 2013; Vaytet et al. 2013). For example, a second core is formed with 4 MJup and 4 R in Masunaga & Inutsuka (2000), while 1.4 MJup and 0.65 R in the ten simulations of Vaytet et al. (2013)10. Therefore, our calculations start with slightly evolved seeds rather than second cores.

Our approach in Sects. 3.4 and 5.1.4 was to fix the initial seed mass and explore different values of the initial seed radius, following HOK11. This choice was essentially motivated by the fact that the convergence of evolution models with MESA is more difficult for very low seed masses. However, in light of the fact that, as pointed out by BVC12, the low value of the initial seed mass has important consequences for the subsequent stellar evolution, we must show that the range of initial conditions yield evolutionary tracks that are equivalent to those that we would obtain with small initial seed masses.

Figure D.1 illustrates seven evolutionary models starting from different seed masses: 0.01 M (I–II), 0.004 M (III–IV), and 10-3M (V–VII), for various initial radii Rini and heat injection efficiencies ξ. Case (I) corresponds to the fiducial initial condition used in the present paper and case (II) is used in Sect. 5.1.4. Case (III) corresponds to seed conditions obtained by Masunaga & Inutsuka (2000). Case (IV) corresponds to a seed with the same specific entropy as our fiducial seed (I) but a mass of only 0.004 M (see Eq. (B.1)). Case (V) corresponds to the largest specific entropy for which we could converge an evolution calculation for a 10-3M seed. The radius of a 10-3M seed with the same entropy as case (I) would be 3.2 R, but unfortunately the corresponding calculation failed to converge. Finally, cases (VI) and (VII) are obtained by arbitrarily changing the initial radii for 10-3M seeds to span a range of specific entropies. Interestingly, for stars evolving from the initial conditions (II), (V), (VI), or (VII), their radius is so small that they reach high enough central temperatures to ignite hydrogen burning when their mass reaches ≃ 0.6 M.

thumbnail Fig. D.1

Evolution of radius as a function of mass for accreting stars with = 10-5M/ yr for different initial conditions and different values of ξ. The crosses (with corresponding I to VII labels) correspond to the initial conditions chosen to span a range of specific entropies. The different initial radii are 1.5 R (case I), 0.25 R (II), 4.0 R (III), 2.0 R (IV), 2.5 R (V), 1.8 R (VI), and 1.1 R (VII). Solid, dashed, and dotted lines indicate the evolution with ξ = 0,0.1 and 0.5, respectively. The evolution for cases (I) and (II) correspond to those adopted in the present manuscript and are highlighted with thicker lines.

We can see in Fig. D.1 that the range of radii obtained for models with different initial seed masses and different values of ξ is the same as that obtained with the our fiducial seed mass and various values of the initial radius. Furthermore, we see that the evolution is largely controlled by the value of the ξ parameter: For ξ = 0.5, the evolutionary tracks converge independent of the choice of the initial seed properties above about 0.1 M. This is also the case for the ξ = 0.1 case, although differences in radii on the order of ~ 10% remain even after accretion is completed. The cold accretion case is the one for which differences in initial conditions persist the longest and can still be on the order of ~ 40% at the end of the accretion phase. Our conclusion that most stars would have been formed with ξ ≳ 0.1 is not affected by the uncertainties on the initial conditions.

In order to explain the small sizes of young cool stars, HOK11 invoked the need for small seed radius (0.2 R at 0.01 M) while BVC12 invoked the need for a small seed mass. Our calculations in Fig. D.1 show that these conditions are equivalent.

All Tables

Table 1

Input parameters.

Table A.1

Results of the χ2 test to reproduce solar values.

All Figures

thumbnail Fig. 1

Evolution of radius (solid line) and mass (dashed) in the cold accretion limit (ξ = 0) and under the fiducial conditions listed in Table 1. The five main evolution phases are labeled as follows: (I) the contraction phase; (II) deuterium-burning phase; (III) second contraction phase until the accretion is completed; (IV) swelling phase; and (V) main sequence from ~ 2 × 107 yr.

In the text
thumbnail Fig. 2

Radius evolution (top panel) and evolutionary tracks in the H-R diagram (bottom) comparing the cold accretion (solid line) with the classical 1 M evolution (dashed). The filled square indicates the initial location, the triangle indicates the stellar mass at 0.1 M, and the circle indicates the end of accretion. The non-accreting classical evolution, whose initial radius is assumed to be 4.92 R as shown by the asterisks, starts at 0.99 × 105 yr to compare the evolution. At 105 yr, the radii are about one order of magnitude different. In the classical model the star evolves along the Hayashi and Henyey tracks, whereas in the cold accretion model the evolution is roughly horizontal in the H-R diagram. Although the deuterium burning slightly increases the luminosity where log Teff ≃ 3.5, it is much smaller than the Henyey track of the non-accreting 1 M star.

In the text
thumbnail Fig. 3

Radius evolution (top panel) and evolutionary tracks (bottom) with varying deuterium contents ranging from XD = 40 (red solid line), 30 (green dashed), 20 (blue dotted), 10 (pink dot-dashed), and 0 (blue double dot-dashed) ppm. The points represent the end of the accretion phase and the black solid lines represent the classical PMS evolutionary tracks.

In the text
thumbnail Fig. 4

Evolutionary tracks for different initial radius Rini and different initial deuterium content XD assuming cold accretion (ξ = 0). The solid, dashed, dotted and dot-dashed lines represent the cases Rini = 0.25,0.5,1.5, and 3 R, respectively. The colors indicate deuterium content, XD = 20 ppm (blue) and XD = 35 ppm (red). The points represent the observed PMS stars in clusters (see Fig. 8). In all cases, stellar masses reach about 0.5 M at 3500 K.

In the text
thumbnail Fig. 5

Radius evolution (top panel) and evolutionary tracks (bottom) with varying final masses ranging from 0.05 M (solid line), 0.1 M (dashed), 0.3 M (dotted), 1 M (dot-dashed), and 1.5 M (double dot-dashed). The accretion is cold (ξ = 0) and steady (10-5M/ yr). The filled circles indicate the ages when the accretion is completed and the pluses, crosses, and asterisks represent 1, 5, and 10 million years, respectively. The calculation stops at 1010 yr or when stars start leaving the MS (only in the case of 1.5 M). We note that the evolutionary track of 0.05 M is cut short.

In the text
thumbnail Fig. 6

Cold-accretion (ξ = 0) evolution of the stellar radius with time for the studies of HOK11 (dashed green) and BC10 (dashed red) and for our results with steady accretion and XD = 35 ppm (plain green), XD = 20 ppm (dotted red). A model with episodic accretion (as described in Sect. 3.2) and XD = 20 ppm is also shown (plain red).

In the text
thumbnail Fig. 7

Top panel: radius evolution for various heat injection efficiencies, ξ = 0.5 (solid line), 0.1 (dashed), 0.05 (dotted), 0.01 (dot-dashed), and 0 (double dot-dashed). The classical PMS evolution, which is indicated by the thin solid line and starts at the asterisks, has a substantial overlap with the lines of ξ = 0.5. Bottom panel: evolutionary tracks with different ξ ranging 0.5 (red line), 0.1 (green), 0.05 (blue), 0.01 (magenta) and 0 (cyan). The solid, dotted, and dashed lines represent the cases with uniform heat injection, mke = 1, and mke = 0.1, respectively. In both panels, we adopt XD = 20 ppm and Rini = 1.5 R.

In the text
thumbnail Fig. 8

Comparison of 0.3 (solid line), 1 (dashed), 3 (dotted), and 10 (dot-dashed) million-year isochrones for the cases ξ = 0.5, 0.1, 0.05, and 0 from top to bottom in the H-R diagram with observed PMS stars: ONC (red, small diamonds), ρ Oph (red squares), σ Ori (green circles), Tau-Aur (blue triangles), and Tau and Cha I (blue, inverted triangles). Each isochrone is comprised of evolutionary tracks of 0.01, 0.25, 0.05, 0.1, 0.3, 0.5, 0.8, and 1 M stars denoted by pluses, crosses, asterisks, open squares, open circles, open triangles, open inverted triangles, and open tilted squares, respectively. The initial radius is 1.5 R and XD = 20 ppm. The isochrones with varying ξ cover the wide range of the region where the observed stars exist in the high-temperature regime.

In the text
thumbnail Fig. 9

Comparison of 1 Myr isochrones for the cases XD = 35 (red lines) and 20 ppm (blue) ranging ξ from 0 (dotted) to 0.5 (solid), respectively. The initial radius is 1.5 R. The points for observation and isochrones are the same as in Fig. 8. The isochrone with (XD) = (35 ppm,0) is almost the same as that with (20 ppm, 0.1), implying that the different deuterium content has a large impact on the PMS evolution. The thin solid line is a 1 Myr non-accreting isochrone (Baraffe et al. 1998).

In the text
thumbnail Fig. 10

Comparison of isochrones at 0.3, 1, 3, and 10 Myr from top to bottom. Rini is either 3 R (gray lines) or 0.25 R (black). The solid, dashed, and dotted lines represent isochrones with ξ = 0.5,0.1, and 0, respectively. Here XD = 20 ppm. The points are the same as in Fig. 8.

In the text
thumbnail Fig. 11

Evolutionary tracks whose final masses are 0.5, 0.8, and 1 M for the cases of ξ = 0 (green lines), 0.05 (blue), and 0.1 (magenta dot-dashed). The initial radius is 1.5 R and XD = 20 ppm (solid lines) and 35 ppm (dashed). Some are merged to the Hayashi track of the non-accreting evolutionary tracks (black dotted), and others are not. Thus, the mass determination of stars on the Hayashi phase is not affected by the variation in accretion heat. As already discussed, the age determination is affected by ξ. Pluses, crosses, and asterisks denote 1, 5, and 10 million years, respectively. The filled and open circles represent the primary and secondary of the eclipsing binaries CoRoT 223992193, respectively (Stassun et al. 2014). For the cases of XD = 20 ppm, their age is estimated to 1 or 5 million years if ξ = 0.05 or 0.1, respectively.

In the text
thumbnail Fig. B.1

Evolution of the intrinsic luminosity (L, solid line), the energy production rate by thermonuclear reaction (Lnuc, dashed) and the seventh part of the gravitational energy of the accreting material (17Lacc\hbox{$\frac{1}{7}L\sub{acc}$}, dotted). The shaded regions and labels from (I) to (V) are the same as in Fig. 1. During phase (II), Lnuc exceeds 17Lacc\hbox{$\frac{1}{7}L\sub{acc}$} and the star expands. The flat Lnuc during phase (III) results from the steady and instantaneous fusion of the newly accreted deuterium.

In the text
thumbnail Fig. B.2

Profiles of luminosity (top panel), entropy (middle), and radius (bottom) in the stellar interior at 0.1 (red solid line; just after the accretion is completed), 1.0 (green dashed), 2.9 (blue dotted), and 5.1 (magenta dot-dashed) million years. As the luminosity propagates toward the surface, the entropy is redistributed. Crosses represent the base of the surface convective zone; kB and mamu are the Boltzmann constant and the atomic mass unit, respectively.

In the text
thumbnail Fig. D.1

Evolution of radius as a function of mass for accreting stars with = 10-5M/ yr for different initial conditions and different values of ξ. The crosses (with corresponding I to VII labels) correspond to the initial conditions chosen to span a range of specific entropies. The different initial radii are 1.5 R (case I), 0.25 R (II), 4.0 R (III), 2.0 R (IV), 2.5 R (V), 1.8 R (VI), and 1.1 R (VII). Solid, dashed, and dotted lines indicate the evolution with ξ = 0,0.1 and 0.5, respectively. The evolution for cases (I) and (II) correspond to those adopted in the present manuscript and are highlighted with thicker lines.

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.