Revisiting the pre-main-sequence evolution of stars I. Importance of accretion efficiency and deuterium abundance
Masanobu Kunitomo Tristan Guillot Taku Takeuchi Shigeru Ida
 


   

   I. Importance of accretion efficiency and deuterium abundance
   
The isochrones are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via                                  http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/599/A49http://cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/599/A49
   
   


   

   Department of Physics, 
                        Nagoya University, Furo-cho, Chikusa-ku, Nagoya, Aichi 464-8602, Japan

                        kunitomo@nagoya-u.jp
                Université de Nice-Sophia Antipolis,
             Observatoire de la Côte d'Azur,
             CNRS UMR 7293,
             06304 Nice CEDEX 04, France
                Department of Earth and Planetary Sciences, 
                         Tokyo Institute of Technology,
                         2-12-1 Ookayama, Meguro-ku, Tokyo 152-8551, Japan
         Earth-Life Science Institute,
                         Tokyo Institute of Technology,
                         2-12-1 Ookayama, Meguro-ku, Tokyo 152-8551, Japan

   



 
  
  
   
   Protostars grow from the first formation of a small seed and subsequent accretion of material.
   Recent theoretical work has shown that the pre-main-sequence (PMS) evolution of stars is much more complex than previously envisioned.
        Instead of the traditional steady, one-dimensional solution,       
   accretion may be episodic and not necessarily symmetrical, thereby affecting the energy deposited inside the star and its interior structure.
   
  
   Given this new framework, we want to understand what controls the evolution of accreting stars.
   
  
   We use the MESA stellar evolution code with various sets of conditions. In particular, we account for the (unknown) efficiency of accretion in burying gravitational energy into the protostar through a parameter, ξ, and we vary the amount of deuterium present.
  
   We confirm the findings of previous works that, in terms of evolutionary tracks on the Hertzsprung-Russell (H-R) diagram, the evolution changes significantly with the amount of energy that is lost during accretion. 
   We find that deuterium burning also regulates the PMS evolution. 
   In the low-entropy accretion scenario, the evolutionary tracks in the H-R diagram are significantly different from the classical tracks and are sensitive to the deuterium content. A comparison of theoretical evolutionary tracks and observations allows us to exclude some cold accretion models (ξ∼ 0) with low deuterium abundances.
   
  
   We confirm that the luminosity spread seen in clusters can be explained by models with a somewhat inefficient injection of accretion heat. The resulting evolutionary tracks then become sensitive to the accretion heat efficiency, initial core entropy, and deuterium content. In this context, we predict that clusters with a higher D/H ratio should have less scatter in luminosity than clusters with a smaller D/H. 
        Future work on this issue should include radiation-hydrodynamic simulations to determine the efficiency of accretion heating and further observations to investigate the deuterium content in star-forming regions.
   

   

   Revisiting the pre-main-sequence evolution of stars
    Masanobu Kunitomo<ref>,
          Tristan Guillot<ref>,
          Taku Takeuchi,<ref>
          Present affiliation: Sanoh Industrial Co., Ltd., Japan 
         Shigeru Ida<ref>
          
    Received 5 February 2016 / Accepted 6 December 2016
================================================================================================================================================================================================







§ INTRODUCTION
 
Since the pioneering work of <cit.>, 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 <cit.>.

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 <cit.>. This has strong consequences for the stellar evolutionary tracks and therefore the inferred physical properties of the star <cit.>. Finally, it controls the level to which planet formation affects stellar surface composition, as suggested from observations <cit.> or theoretical models <cit.>.

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. <ref>, we describe our physical model and computation method for simulating the PMS evolution including accretion. 
In Sect. <ref>, 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. <ref>, we explore the dependence on the entropy of accreting matter.
In Sect. <ref>, we compile the results of Sects. <ref> and <ref> to evaluate observational constraints and observational consequences. Our results are summarized in Sect. <ref>.





§ METHOD





 §.§ 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 <cit.> 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 <cit.> and <cit.>. 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 <cit.>, 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 <cit.> (see Appendix <ref>).
The energy equation is described in Sect. <ref> in detail.

We use the radiative opacity of <cit.> and <cit.> for the temperature range below 10^4.5 K, and the electron conduction opacity of <cit.>.
We adopt the MESA default version of the equation of state; we basically follow the OPAL EOS tables <cit.> and SCVH tables <cit.>.






 §.§ Thermonuclear reactions
 
Pre-main-sequence stars ignite light elements such as deuterium and lithium when the central temperature exceeds ∼ 10^6 and ∼3×10^6 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 ∼ 10^3 times smaller abundance <cit.>.
We use the thermonuclear reaction rates of <cit.>.
Main-sequence (MS) stars ignite hydrogen when temperature exceeds ∼ 10^7 K.






 §.§ 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 <ref>). 
We conducted a χ^2 test and used the following converged input parameters: X_ini=0.70046, Y_ini=0.27948 and Z_ini=0.02006.
We assume ^3He/^4He = 10^-4, similar to the value in the Jovian atmosphere <cit.>.
We use the composition of metals of <cit.>.


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. <ref>.

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.






 §.§ Initial conditions


Stars are formed by the collapse of molecular cloud cores through a number of stages <cit.>. 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 () 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, respectively, following <cit.> and <cit.>. 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 <cit.>.
Therefore initial seeds in the present work correspond to slightly evolved protostars from second cores.
As discussed by <cit.>, 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 <ref> that for ranges of masses between 1 and 4, the range of radii when the mass attains 0.01 M_⊙ is between 0.25 and 1.5. Our fiducial value of 1.5 therefore corresponds to a high initial value of the entropy of the star. In Sect. <ref> we explore the consequences of variations of the initial radius down to 0.25.





 §.§ Mass accretion
 

The collapse of the molecular cloud core yields an accretion rate Ṁ∼ cs^3/G, where cs is the characteristic speed of sound in the cloud <cit.>. For molecular cloud temperatures T∼ 30 K, this implies Ṁ≈ 10^-5  M_⊙/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 <cit.>. 
It can be strongly variable (episodic), as observed in FU Ori.


Our fiducial accretion rate is Ṁ≈ 10^-5  M_⊙/yr but we explore the effect of varying this rate in Sect. <ref>.
We stop mass accretion abruptly when the stellar mass reaches the final mass, Mfin. 





 §.§ 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 <cit.>. 

First, we assume that the specific entropy of the accreting material is the same as that at the stellar surface. 
This follows from <cit.> 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_⋆,

where M_⋆ and R_⋆ are stellar mass and radius, respectively, and energy constraints impose that ξ is such that 0≤ξ≤ 1
[
Our ξ is equivalent to αϵ in <cit.>.
].
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, <cit.> 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.


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 <cit.>, and in which Ladd is distributed uniformly and instantaneously within the entire star. The energy deposited per unit mass is thus 

    εadd^ (uniform)=Ladd/M_⋆.

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)=Ladd/M_⋆ Max[0,2/mke^2( M_r/M_⋆-1+mke) ].

where M_r 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 ∫εadd dM=Ladd. The energy conservation equation thus is written 

    ∂ L/∂ M_r = ε_nuc-( T∂ S/∂ t)_M_r+ε_add,

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  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) 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 <cit.>. 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,  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 <cit.>. 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 <cit.>.
MESA originally treated the thermodynamic state of the accreting material with compressional heating <cit.>. 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 <cit.>.









 §.§ 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 and R_⋆ > 0.7, 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 <cit.> and ATLAS9 <cit.>. 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 <cit.> model atmospheres.
The impact of the switch of the boundary condition on the stellar structure and evolution is negligible.




 §.§ Adopted parameters and comparison with previous studies
 

Table <ref> 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. <ref>–<ref> 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 ( and ) in Sect. <ref>.







§ 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. <ref>).


Figure <ref> shows how the radius changes with time in the cold accretion case and with our fiducial settings (see Table <ref>).
As described in more detail in Appendix <ref>, 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 × 10^7 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. <ref>), we examine the influence of mass accretion rate, deuterium content, final stellar mass, and initial conditions in Sects. <ref>–<ref>.
We compare our results against previous studies in Sect. <ref>.













 §.§ Cold versus classical models

We now compare the radius evolution in the case of cold steady accretion with that in the classical model.
Figure <ref>a 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 <cit.>, we adopt for the classical evolution case an initial radius of 4.92 at 10^5 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.






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. <ref>b.
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 Hertzsprung-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. <ref>.

We point out that the evolutionary track in the cold accretion case in Fig. <ref>b 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 .



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 <ref>.





 §.§ Dependence on mass accretion rate


As described in Sect. <ref>, 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^-6 /yr,
and time variability; i.e., the episodic accretion following the pioneering works of ,  and .
For episodic accretion, we adopt the parameters of  on the basis of hydrodynamic simulations by <cit.>.


We confirm the findings of  and  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 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^-5 /yr from now on.






 §.§ Dependence on deuterium content


The PMS evolution is controlled by deuterium burning when it occurs vigorously (i.e., during phase II in Sect. <ref>).
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 <cit.> who concluded that the radius is only moderately affected by deuterium content.
For example, when stars become 1 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 content[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).]. 

The primordial value at the beginning of the universe is XD,prim=35.8±2.5 ppm <cit.>, 

the indirectly[
Two methods may be used to constrain the protosolar deuterium abundance: using either the deuterium content in the Jovian atmosphere <cit.> or the enhanced ^3He measured in the solar wind and compared to the Jovian atmosphere <cit.>. 
]
estimated value of the protosolar nebula is XD,PSN=28.0±2.8 ppm <cit.>, 

and the present-day value of the local interstellar medium varies widely from XD,ISM=13.7±5.3 ppm <cit.> to at least 32.3 ± 3.4 <cit.>.
Classical studies, such as <cit.>, used a higher value (XD=35 ppm) based on classical observations of local interstellar media <cit.>.
 
Although these values still remain uncertain to some extent and, in particular, XD,ISM is still under debate <cit.>, 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 <ref> 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 1/7Lacc for a longer period of time and then the star becomes larger (see Eq. <ref>), where Lacc≡ GM_⋆Ṁ/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.






In the extreme case XD=0, the star keeps shrinking even after ∼ 10^4 years.
Just before the accretion is completed at 10^5 years, the radius becomes ∼ 0.3 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
<cit.> with spherical accretion. 
This is because the total energy injected by accretion, Eacc=∫ Ladddt, exceeds the total energy generated by deuterium burning, 
ED=qD XD Mfin, where qD=2.63× 10^18  erg g^-1 is the energy released by deuterium fusion per gram of deuterium, when


    ξ≳ 0.05 XD20 ppmR̅Mfin^-1,


where 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. 








 §.§ Dependence on initial conditions


So far we have assumed that our initial stellar seed is characterized by 1.5 and 0.01. As pointed out by , since the entropy of accreting materials is related to the stellar entropy in our prescription (see Sect. <ref>), the initial conditions affect the evolution <cit.>. This is also pointed out by  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 <ref>, the approach is equivalent with the current uncertainties in stellar seed masses leading to values of the radius at 0.01 between 0.25 and 1.5.




In Fig. <ref>, we show the evolutionary tracks of accreting protostars with brown-dwarf masses obtained with Rini from 0.25 to 3 and XD=20 ppm and 35 ppm. We do not show the results for Rini≳3 since these are similar to the Rini=3 case because of their short K-H timescale.
Like  (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 for Rini=1.5 and 0.25, respectively; see Fig. <ref>).








We see that the evolutionary track with Rini=0.25 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. <ref>).
Moreover, the effect of deuterium burning has less impact in this case. 
The star expands in Teff≃2500–2900 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, Lacc[The luminosity in the present paper is always the stellar intrinsic luminosity and does not include the accretion luminosity emitted from the shock surface.].
As described in Appendix <ref>, stars expand when LD>1/7Lacc (Eq. <ref>). Hence the larger Lacc makes the expansion phase (phase II in Sect. <ref>) shorter and the stellar radius and luminosity remain small (see also ).









 §.§ 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 <ref>.
The cold accretion evolution of stars with masses of 0.05, 0.1, 0.3, 1, and 1.5 are shown in Fig. <ref>.
Low-mass stars contract for a long time because of their long K-H timescale. For example, stars with 0.1 and 0.3 enter their MS at 10^9 and 3 × 10^8 years, respectively. During the first million years of evolution, Fig. <ref> shows that the size of stars undergoing cold accretion is not a monotonic function of mass. For example, a 0.1 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 star contracts monotonically and even more after ≃ 1–2 million years when all deuterium has been consumed. 
Higher mass stars (1.0, 1.5) expand just after accretion has stopped because of hydrogen burning. 













 §.§ Comparison with previous studies
 

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


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. <ref>). 
Figure <ref> compares the radius evolution in  and  to our calculations with the same hypotheses. Given the fact that different stellar models are used (including different values of αMLT –see Table <ref>), there is excellent agreement between our calculations and those of  and . 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.
















§ 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. (<ref>). Then, we examine, with the parameter mke (see Eq. <ref>), 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 <ref>).




 §.§ Early evolution with a varying ξ


Figure <ref> 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. 

The evolutionary track of accreting stars corresponds to the “birthline” proposed by <cit.>.
As  pointed out, Fig. <ref>b 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 10^5 yr) in the case of ξ=0.5 is almost the same as obtained by <cit.> for spherical accretion.










Figure <ref>a  shows that before deuterium burning (at 10^4 to 10^5 yrs), 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 <ref>.
The situation corresponds to a case when nuclear reactions and radiative cooling may be neglected (i.e., Lnuc, L_⋆≪ Lacc). The energy equation may thus be shown to yield


    Ṙ/R_⋆ = ( 2-1-ξ/C) Ṁ/M_⋆ 
    ,


where C is a constant as a function of the specific heat ratio and the polytropic index (see Appendix <ref>).
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×10^5 years, 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 <ref> and Eq. <ref>).









 §.§ Importance of the location of heat injection

So far, we assumed a uniform distribution of injected heat as in  (see Eq. <ref>).
However, this assumption may not be valid, especially in stars with a large radiative core.
Moreover, as described in Sect. <ref>, recent radiation-hydrodynamic simulations by <cit.> 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. <ref>).

Figure <ref>b 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. <ref>).
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. (<ref>), 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. (<ref>) is a useful simplification to represent possible evolutionary tracks.







§ 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. <ref>.
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 (,  and ), we choose here to attempt to constrain the values of ξ that are in agreement with the observations of young clusters. 







 §.§ The PMS luminosity spread problem




  §.§.§ Observational constraints

The luminosity spread of PMS stars has been a matter of debate for decades <cit.>.
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 <cit.>, (ii) an intrinsic age spread <cit.> and (iii) physical processes <cit.>.
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 <cit.>.
However, the contribution to the luminosity spread by each uncertainty has been quantitatively estimated <cit.> and various authors claimed that the sum of their contributions is smaller than the observed luminosity spread. Moreover, <cit.> 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, <cit.> 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 ;  see also Sect. <ref>). For stars of lower temperatures, the results depend on the size of the assumed seed radius and differ between the different studies (see Sect. <ref>). 




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. <ref>).







  §.§.§ Isochrones as a function of assumed ξ
 

In Fig. <ref> 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, and Mini=0.01).
We use the following observational data of young stars: ρ Ophiuchus <cit.>, σ Orionis <cit.>, Taurus and Chamaeleon I <cit.>, Taurus-Auriga <cit.>, and Orion nebula cluster <cit.>.








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. <cit.> 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 , , and .
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. <ref>.







  §.§.§ Dependence on deuterium content
 

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

Figure <ref> 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. <ref>, this is because in the case of cold accretion, deuterium burning regulates the PMS evolution.
But Fig. <ref> 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. 
















  §.§.§ Dependence on initial conditions
 

As described in Sect. <ref> 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. <ref>, we show the isochrones obtained with Rini of 0.25 and 3, 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. <ref>) and corresponds to the situation described in Sect. <ref>. 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.







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




  §.§.§ 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 ∼ 1Myr[
We do not derive the age of Tau and Cha I because the data in <cit.> is an assembly of young stars in two star-forming regions.].
However, inaccurate estimations of stellar luminosities, membership issues, and other problems can affect these age estimates. 






We can see in Sect. <ref> and Fig. <ref> that for a majority of stars ξ>0.1 for XD=20 ppm and Rini=1.5 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, 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. <ref>). 








 §.§ The case of CoRoT 223992193


We now consider the case of a specific system, the eclipsing binary CoRoT 223992193 <cit.>. 
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 <ref> 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. 
The two stars have masses between 0.6 and 0.8 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 ξ. 









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.








§ 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). 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 <cit.> and <cit.> 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) 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, 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. 





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.




aa
86
natexlab#1#1

[Allard et al.(2001)Allard, Hauschildt, Alexander,
  Tamanai, & Schweitzer]Allard+01
Allard, F., Hauschildt, P. H., Alexander, D. R., Tamanai, A., &
  Schweitzer, A. 2001, , 556, 357

[Amelin et al.(2002)Amelin, Krot, Hutcheon, &
  Ulyanov]Amelin+02
Amelin, Y., Krot, A. N., Hutcheon, I. D., & Ulyanov, A. A. 2002,
  Science, 297, 1678

[Angulo et al.(1999)Angulo, Arnould, Rayet, Descouvemont,
  Baye, Leclercq-Willain, Coc, Barhoumi, Aguer, Rolfs, Kunz,
  Hammer, Mayer, Paradellis, Kossionides, Chronidou, Spyrou,
  degl'Innocenti, Fiorentini, Ricci, Zavatarelli, Providencia,
  Wolters, Soares, Grama, Rahighi, Shotter, & Lamehi
  Rachti]Angulo+99
Angulo, C., Arnould, M., Rayet, M., et al. 1999, Nuclear Physics A,
  656, 3

[Asplund et al.(2009)Asplund, Grevesse, Sauval, &
  Scott]Asplund+09
Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, , 47,
  481

[Bahcall et al.(2005)Bahcall, Basu, Pinsonneault, &
  Serenelli]Bahcall+05
Bahcall, J. N., Basu, S., Pinsonneault, M., & Serenelli, A. M. 2005,
  , 618, 1049

[Baraffe & Chabrier(2010)]BC10
Baraffe, I. & Chabrier, G. 2010, , 521, A44

[Baraffe et al.(1998)Baraffe, Chabrier, Allard, &
  Hauschildt]Baraffe+98
Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, ,
  337, 403

[Baraffe et al.(2009)Baraffe, Chabrier, & Gallardo]BCG09
Baraffe, I., Chabrier, G., & Gallardo, J. 2009, , 702, L27

[Baraffe et al.(2012)Baraffe, Vorobyov, & Chabrier]BVC12
Baraffe, I., Vorobyov, E., & Chabrier, G. 2012, , 756, 118

[Basu & Antia(2004)]Basu+Antia04
Basu, S. & Antia, H. M. 2004, , 606, L85

[Burningham et al.(2005)Burningham, Naylor, Littlefair, &
  Jeffries]Burningham+05
Burningham, B., Naylor, T., Littlefair, S. P., & Jeffries, R. D. 2005,
  , 363, 1389

[Cassisi et al.(2007)Cassisi, Potekhin, Pietrinferni,
  Catelan, & Salaris]Cassisi+07
Cassisi, S., Potekhin, A. Y., Pietrinferni, A., Catelan, M., &
  Salaris, M. 2007, , 661, 1094

[Castelli & Kurucz(2004)]Castelli+Kurucz04
Castelli, F. & Kurucz, R. L. 2004, ArXiv Astrophysics e-prints

[Chabrier et al.(2007)Chabrier, Gallardo, &
  Baraffe]Chabrier+07
Chabrier, G., Gallardo, J., & Baraffe, I. 2007, , 472, L17

[Chambers(2010)]Chambers10
Chambers, J. E. 2010, , 724, 92

[Chandrasekhar(1967)]Chandrasekhar67
Chandrasekhar, S. 1967, An introduction to the study of stellar structure

[Cox & Giuli(1968)]Cox+Giuli68
Cox, J. & Giuli, R. 1968, Gordon and Breach, New York, 401

[Da Rio et al.(2010)Da Rio, Robberto, Soderblom, Panagia,
  Hillenbrand, Palla, & Stassun]Da-Rio+10
Da Rio, N., Robberto, M., Soderblom, D. R., et al. 2010, , 722,
  1092

[Dunham & Vorobyov(2012)]Dunham+Vorobyov12
Dunham, M. M. & Vorobyov, E. I. 2012, , 747, 52

[Ferguson et al.(2005)Ferguson, Alexander, Allard, Barman,
  Bodnarik, Hauschildt, Heffner-Wong, & Tamanai]Ferguson+05
Ferguson, J. W., Alexander, D. R., Allard, F., et al. 2005, , 623,
  585

[Gatti et al.(2008)Gatti, Natta, Randich, Testi, &
  Sacco]Gatti+08
Gatti, T., Natta, A., Randich, S., Testi, L., & Sacco, G. 2008,
  , 481, 423

[Gatti et al.(2006)Gatti, Testi, Natta, Randich, &
  Muzerolle]Gatti+06
Gatti, T., Testi, L., Natta, A., Randich, S., & Muzerolle, J. 2006,
  , 460, 547

[Geroux et al.(2016)Geroux, Baraffe, Viallet, Goffrey,
  Pratt, Constantino, Folini, Popov, & Walder]Geroux+16
Geroux, C., Baraffe, I., Viallet, M., et al. 2016, , 588, A85

[Gillen et al.(2014)Gillen, Aigrain, McQuillan, Bouvier,
  Hodgkin, Alencar, Terquem, Southworth, Gibson, Cody, Lendl,
  Morales-Calderón, Favata, Stauffer, & Micela]Gillen+14
Gillen, E., Aigrain, S., McQuillan, A., et al. 2014, , 562, A50

[Grevesse & Sauval(1998)]GS98
Grevesse, N. & Sauval, A. J. 1998, , 85, 161

[Guillot(1999)]Guillot99
Guillot, T. 1999, , 47, 1183

[Guillot et al.(2014)Guillot, Ida, & Ormel]Guillot+14
Guillot, T., Ida, S., & Ormel, C. W. 2014, , 572, A72

[Hartmann(2001)]Hartmann01
Hartmann, L. 2001, , 121, 1030

[Hartmann et al.(1997)Hartmann, Cassen, &
  Kenyon]Hartmann+97
Hartmann, L., Cassen, P., & Kenyon, S. J. 1997, , 475, 770

[Hauschildt et al.(1999a)Hauschildt, Allard, &
  Baron]Hauschildt+99a
Hauschildt, P. H., Allard, F., & Baron, E. 1999a, ,
  512, 377

[Hauschildt et al.(1999b)Hauschildt, Allard,
  Ferguson, Baron, & Alexander]Hauschildt+99b
Hauschildt, P. H., Allard, F., Ferguson, J., Baron, E., & Alexander,
  D. R. 1999b, , 525, 871

[Hayashi(1961)]Hayashi61
Hayashi, C. 1961, , 13, 450

[Heber et al.(2008)Heber, Baur, Bochsler, Burnett,
  Reisenfeld, Wieler, & Wiens]Heber+08
Heber, V. S., Baur, H., Bochsler, P., et al. 2008, in Lunar and
  Planetary Inst. Technical Report, Vol. 39, Lunar and Planetary Science
  Conference, 1779

[Hébrard et al.(2005)Hébrard, Tripp, Chayer,
  Friedman, Dupuis, Sonnentrucker, Williger, & Moos]Hebrard+05
Hébrard, G., Tripp, T. M., Chayer, P., et al. 2005, , 635, 1136

[Henyey et al.(1965)Henyey, Vardya, &
  Bodenheimer]Henyey+65
Henyey, L., Vardya, M. S., & Bodenheimer, P. 1965, , 142, 841

[Herwig(2000)]Herwig00
Herwig, F. 2000, , 360, 952

[Hillenbrand(2009)]Hillenbrand09
Hillenbrand, L. A. 2009, in IAU Symposium, Vol. 258, IAU Symposium, ed. E. E.
  Mamajek, D. R. Soderblom, & R. F. G. Wyse, 81–94

[Hosokawa et al.(2011)Hosokawa, Offner, &
  Krumholz]Hosokawa+11
Hosokawa, T., Offner, S. S. R., & Krumholz, M. R. 2011, , 738, 140

[Hosokawa & Omukai(2009)]Hosokawa+Omukai09
Hosokawa, T. & Omukai, K. 2009, , 691, 823

[Hueso & Guillot(2005)]Hueso+Guillot05
Hueso, R. & Guillot, T. 2005, , 442, 703

[Inutsuka(2012)]Inutsuka12
Inutsuka, S.-i. 2012, Progress of Theoretical and Experimental Physics, 2012,
  010000

[Inutsuka et al.(2015)Inutsuka, Inoue, Iwasaki, &
  Hosokawa]Inutsuka+15
Inutsuka, S.-i., Inoue, T., Iwasaki, K., & Hosokawa, T. 2015, ,
  580, A49

[Inutsuka et al.(2010)Inutsuka, Machida, &
  Matsumoto]Inutsuka+10
Inutsuka, S.-i., Machida, M. N., & Matsumoto, T. 2010, , 718, L58

[Jeffries(2007)]Jeffries07
Jeffries, R. D. 2007, , 381, 1169

[Jeffries(2012)]Jeffries12
Jeffries, R. D. 2012, Are There Age Spreads in Star Forming Regions?, ed.
  A. Moitinho & J. Alves, 163

[Kenyon & Hartmann(1995)]Kenyon+Hartmann95
Kenyon, S. J. & Hartmann, L. 1995, , 101, 117

[Larson(1969)]Larson69
Larson, R. B. 1969, , 145, 271

[Lellouch et al.(2001)Lellouch, Bézard, Fouchet,
  Feuchtgruber, Encrenaz, & de Graauw]Lellouch+01
Lellouch, E., Bézard, B., Fouchet, T., et al. 2001, , 370, 610

[Linsky et al.(2006)Linsky, Draine, Moos, Jenkins, Wood,
  Oliveira, Blair, Friedman, Gry, Knauth, Kruk, Lacour, Lehner,
  Redfield, Shull, Sonneborn, & Williger]Linsky+06
Linsky, J. L., Draine, B. T., Moos, H. W., et al. 2006, , 647, 1106

[Mahaffy et al.(1998)Mahaffy, Donahue, Atreya, Owen, &
  Niemann]Mahaffy+98
Mahaffy, P. R., Donahue, T. M., Atreya, S. K., Owen, T. C., &
  Niemann, H. B. 1998, , 84, 251

[Martin et al.(2012)Martin, Lubow, Livio, &
  Pringle]Martin+12
Martin, R. G., Lubow, S. H., Livio, M., & Pringle, J. E. 2012, ,
  423, 2718

[Masunaga & Inutsuka(2000)]Masunaga+Inutsuka00
Masunaga, H. & Inutsuka, S.-i. 2000, , 531, 350

[Meléndez et al.(2009)Meléndez, Asplund, Gustafsson,
  & Yong]Melendez+09
Meléndez, J., Asplund, M., Gustafsson, B., & Yong, D. 2009, ,
  704, L66

[Mercer-Smith et al.(1984)Mercer-Smith, Cameron, &
  Epstein]Mercer-Smith+84
Mercer-Smith, J. A., Cameron, A. G. W., & Epstein, R. I. 1984, ,
  279, 363

[Muzerolle et al.(2005)Muzerolle, Luhman, Briceño,
  Hartmann, & Calvet]Muzerolle+05
Muzerolle, J., Luhman, K. L., Briceño, C., Hartmann, L., &
  Calvet, N. 2005, , 625, 906

[Nelder & Mead(1965)]Nelder+Mead65
Nelder, J. A. & Mead, R. 1965, The computer journal, 7, 308

[Palla & Stahler(1992)]Palla+Stahler92
Palla, F. & Stahler, S. W. 1992, , 392, 667

[Palla & Stahler(2000)]Palla+Stahler00
Palla, F. & Stahler, S. W. 2000, , 540, 255

[Paxton et al.(2011)Paxton, Bildsten, Dotter, Herwig,
  Lesaffre, & Timmes]Paxton+11
Paxton, B., Bildsten, L., Dotter, A., et al. 2011, , 192, 3

[Paxton et al.(2013)Paxton, Cantiello, Arras, Bildsten,
  Brown, Dotter, Mankovich, Montgomery, Stello, Timmes, &
  Townsend]Paxton+13
Paxton, B., Cantiello, M., Arras, P., et al. 2013, , 208, 4

[Paxton et al.(2015)Paxton, Marchant, Schwab, Bauer,
  Bildsten, Cantiello, Dessart, Farmer, Hu, Langer, Townsend,
  Townsley, & Timmes]Paxton+15
Paxton, B., Marchant, P., Schwab, J., et al. 2015, , 220, 15

[Prantzos(2007)]Prantzos+07
Prantzos, N. 2007, , 130, 27

[Ramírez et al.(2009)Ramírez, Meléndez, &
  Asplund]Ramirez+09
Ramírez, I., Meléndez, J., & Asplund, M. 2009, , 508, L17

[Ramírez et al.(2011)Ramírez, Meléndez,
  Cornejo, Roederer, & Fish]Ramirez+11
Ramírez, I., Meléndez, J., Cornejo, D., Roederer, I. U., &
  Fish, J. R. 2011, , 740, 76

[Reggiani et al.(2011)Reggiani, Robberto, Da Rio, Meyer,
  Soderblom, & Ricci]Reggiani+11
Reggiani, M., Robberto, M., Da Rio, N., et al. 2011, , 534, A83

[Rogers & Nayfonov(2002)]Rogers+Nayfonov02
Rogers, F. J. & Nayfonov, A. 2002, , 576, 1064

[Saumon et al.(1995)Saumon, Chabrier, & van
  Horn]Saumon+95
Saumon, D., Chabrier, G., & van Horn, H. M. 1995, , 99, 713

[Seaton(2005)]Seaton05
Seaton, M. J. 2005, , 362, L1

[Shu(1977)]Shu77
Shu, F. H. 1977, , 214, 488

[Soderblom et al.(2014)Soderblom, Hillenbrand, Jeffries,
  Mamajek, & Naylor]Soderblom+14
Soderblom, D. R., Hillenbrand, L. A., Jeffries, R. D., Mamajek, E. E.,
  & Naylor, T. 2014, Protostars and Planets VI, 219

[Stahler(1983)]Stahler+83
Stahler, S. W. 1983, , 274, 822

[Stahler(1988)]Stahler88
Stahler, S. W. 1988, , 332, 804

[Stahler & Palla(2005)]Stahler+Palla05
Stahler, S. W. & Palla, F. 2005, The Formation of Stars

[Stahler et al.(1986)Stahler, Palla, &
  Salpeter]Stahler+86
Stahler, S. W., Palla, F., & Salpeter, E. E. 1986, , 302, 590

[Stahler et al.(1980)Stahler, Shu, & Taam]SST80I
Stahler, S. W., Shu, F. H., & Taam, R. E. 1980, , 241, 637

[Stamatellos et al.(2007)Stamatellos, Whitworth, Bisbas, &
  Goodwin]Stamatellos+07
Stamatellos, D., Whitworth, A. P., Bisbas, T., & Goodwin, S. 2007,
  , 475, 37

[Stassun et al.(2014)Stassun, Feiden, &
  Torres]Stassun+14
Stassun, K. G., Feiden, G. A., & Torres, G. 2014, , 60, 1

[Steigman(2006)]Steigman06
Steigman, G. 2006, International Journal of Modern Physics E, 15, 1

[Sugimoto & Nomoto(1975)]Sugimoto+Nomoto75
Sugimoto, D. & Nomoto, K. 1975, , 27, 197

[Tomida et al.(2013)Tomida, Tomisaka, Matsumoto, Hori,
  Okuzumi, Machida, & Saigo]Tomida+13
Tomida, K., Tomisaka, K., Matsumoto, T., et al. 2013, , 763, 6

[Townsley & Bildsten(2004)]Townsley+Bildsten04
Townsley, D. M. & Bildsten, L. 2004, , 600, 390

[Vaytet et al.(2013)Vaytet, Chabrier, Audit, Commerçon, Masson, Ferguson, & Delahaye]Vaytet+13
Vaytet, N., Chabrier, G., Audit, E., et al. 2013, , 557, A90

[Vidal-Madjar & Gry(1984)]Vidal-Madjar+Gry84
Vidal-Madjar, A. & Gry, C. 1984, , 138, 285

[Vorobyov & Basu(2005)]VB05
Vorobyov, E. I. & Basu, S. 2005, , 633, L137

[Vorobyov & Basu(2010)]Vorobyov+Basu10
Vorobyov, E. I. & Basu, S. 2010, , 719, 1896

[Winkler & Newman(1980)]Winkler+Newman80
Winkler, K.-H. A. & Newman, M. J. 1980, , 236, 201




 




§ Χ^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 <cit.>.
Although the solar metallicity remains a matter of debate[
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% <cit.>.]
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” <cit.>. 
The input parameters are the initial composition (X_ini, Y_ini, and Z_ini; see Sect. <ref>), the mixing-length parameter, α_MLT, and the overshoot mixing parameter, f_ov <cit.>. 


We calculated the evolution of 1 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 <cit.>.
The results of the χ^2 test listed in Table <ref> are used in this paper.



   











§ 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. <ref>.
The evolution can be split into five phases as shown in Fig. <ref>. 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

    R_⋆∝ M_⋆^-1/3 .

Indeed, we show in Appendix <ref> 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, τacc≡ M_⋆/Ṁ, where Etot is the total energy of the star and L_⋆ is the stellar intrinsic luminosity [
The total luminosity of protostars is the sum of the intrinsic luminosity, L_⋆, and the radiation from the accretion shock front (i.e., Lacc - Ladd).
].
This is equivalent to Lacc≫ L_⋆ as shown in Fig. <ref>.
The accretion luminosity is defined as 


    Lacc   ≡    GM_⋆Ṁ/R_⋆ ,

which is for example 31.3 in the case of M_⋆=0.1, R_⋆=1, and Ṁ=10^-5 /yr.









II. The deuterium-burning phase:

After the central temperature exceeds ∼ 10^6 K, deuterium fusion affects the evolution (see Sect. <ref>). In the current settings, it happens at t = 7 × 10^3 yr and M_⋆≃ 0.08.
As described in Sect. <ref>, the energy production rate of deuterium burning has a strong temperature dependence, that is, εnuc∝ (T/10^6 K)^11.8. 
This strong temperature sensitivity 
is responsible for a rapid expansion of the star through the so-called “thermostat effect” <cit.>. 


Indeed, after the ignition of deuterium fusion, the central temperature is maintained constant at about 10^6 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 <ref>, Tc∝ Pc/ρc∝ M_⋆/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 <ref>, we derive Eq. (<ref>), 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 <cit.>.
By neglecting the intrinsic luminosity (since the Kelvin-Helmholtz timescale is long compared to relevant timescales in this phase),


    Ṙ/R_⋆ = -1/3Ṁ/M_⋆+7/3R_⋆ Lnuc/GM_⋆^2 .


Therefore, the condition for the expansion is given by


    Lnuc > 1/7Lacc .


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 × 10^4 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.65 L_⊙Ṁ10^-5 /yrXD2.0×10^-5 .

This may be compared to the accretion luminosity Lacc defined in Eq. (<ref>), but for a mass that is now M_⋆∼ 0.5. We thus obtain LD,max/Lacc∼ 0.05 for standard values of the parameters. According to Eq. (<ref>), 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 Lnuc≃ LD,max until the accretion ceases at 10^5 yr (see Fig. <ref>).




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 ∼ 10^7–10^8 years. 
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 ∼ 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. (<ref>).

After several million years, the star expands. 
A “luminosity wave” <cit.> 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 <cit.>.
These profiles are shown in Fig. <ref>.








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.






§ ANALYTICAL RELATIONS
 
Here we derive the mass-radius relationship of Eq. (<ref>) in two ways [In addition, another derivation is possible using a polytropic analysis <cit.>.].
First, we use the characteristic density and pressure and the entropy following <cit.>.
The characteristic density and pressure of the star in the hydrostatic equilibrium are given by

    ρ̃   ∝ M_⋆/R_⋆^3 ,    P̃   ∝ M_⋆^2/R_⋆^4  .

The entropy is given by S=c_Vln( P/ρ^γad) +S_0, 
where c_V is the specific heat at constant volume, γad =c_P/c_V the specific heat ratio, and S_0 is the constant. Substituting Eq. (<ref>) into this equation, we obtain

    R_⋆=M_⋆^-2-γad/3γad-4exp( 1/3γad-4S-S_0/c_V) .

In the case of the monatomic ideal gas (γad=5/3), we obtain

    R_⋆= M_⋆^-1/3exp[ 2/3μ/ℛ (S-S_0)   ] ,

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, R_⋆∝ M_⋆^-1/3 (Eq. <ref>).


Secondly, Eq. (<ref>) can be derived from the energy equation, which also gives Eq. (<ref>) following <cit.>.
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:

    dEtot/dt =-L_⋆ +Lnuc -GM_⋆Ṁ/R_⋆+Ladd .

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. (<ref>) here we assume Ladd=ξ GM_⋆Ṁ/R_⋆ . 
From the virial theorem,


    Etot=(4-3γad)Eint=3γad-4/3(γad-1)Eg .


If we assume the polytropic relation that P(ρ)=Kρ^1+1/n, 


    Eg=-3/5-nGM_⋆^2/R_⋆ ,


where n is the polytropic index. 
If we define C≡(3γad-4)/(γad-1)(5-n), Etot=-CGM_⋆^2/R_⋆. 
Thus, the total energy evolution in Eq. (<ref>) is transformed as


    Ṙ/R_⋆ = ( 2-1-ξ/C) Ṁ/M_⋆ -RL_⋆/CGM^2+RLnuc/CGM^2 .



The second term of RHS of Eq. (<ref>) 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 τacc≡ M/Ṁ [This is equivalent to the condition L_⋆≪ Lacc as shown in Fig. <ref>.].
Therefore, the second term of RHS of Eq. (<ref>) 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. (<ref>).
Moreover if Lnuc≪ Lacc, it is transferred as Eq. (<ref>).








§ DEPENDENCE ON THE INITIAL STELLAR SEED MASS
 


In this paper we chose 10^-2 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 
<cit.>.
For example, a second core is formed with 4 and 4 in <cit.>, while 1.4 and 0.65 in the ten simulations of <cit.> [<cit.> indicate that the difference probably results from different opacity tables.].
Therefore, our calculations start with slightly evolved seeds rather than second cores.

Our approach in Sects. <ref> and <ref> was to fix the initial seed mass and explore different values of the initial seed radius, following . 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 , 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 <ref> illustrates seven evolutionary models starting from different seed masses: 0.01 (I–II),  0.004 (III–IV),  and 10^-3 (V–VII), for various initial radii R_ ini and heat injection efficiencies ξ. 
Case (I) corresponds to the fiducial initial condition used in the present paper and case (II) is used in Sect. <ref>.
Case (III) corresponds to seed conditions obtained by <cit.>. Case (IV) corresponds to a seed with the same specific entropy as our fiducial seed (I) but a mass of only 0.004 (see Eq. <ref>).
Case (V) corresponds to the largest specific entropy for which we could converge an evolution calculation for a 10^-3 seed. 
The radius of a 10^-3 seed with the same entropy as case (I) would be 3.2, but unfortunately the corresponding calculation failed to converge. Finally, cases (VI) and (VII) are obtained by arbitrarily changing the initial radii for 10^-3 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.







We can see in  Fig. <ref> 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. 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,  invoked the need for small seed radius (0.2 at 0.01) while  invoked the need for a small seed mass. 
Our calculations in Fig. <ref> show that these conditions are equivalent.








