© ESO, 2016

1. Introduction

Recent analysis of the effect of accretion on the structure and evolution of young low-mass stars and brown dwarfs suggests that early accretion history may impact the properties of objects (luminosity, radius and effective temperature) in star-formation regions and young clusters even after several Myr when accretion process has ceased (Baraffe et al. 2009, 2012; Baraffe & Chabrier 2010; Hosokawa et al. 2011). This scenario can provide an explanation for the luminosity spread observed in young clusters, without requiring an age spread (Baraffe et al. 2009, 2012). More and more observational evidence suggests that low-mass stars accrete most of their material episodically, through bursts of various intensities. Numerical simulations and analytical studies confirm the existence of such accretion bursts, with different possible driving mechanisms invoked (see Audard et al. 2014, for a comprehensive review). At the very early stages of pre-stellar core collapse, the accretion bursts, as they release a significant amount of energy, can have a profound effect on the thermal and chemical evolution of the protostellar disk and the envelope. They heat up the disk and the inner envelope, raising the gas temperatures by a factor of two or more (Vorobyov et al. 2014) and evaporating CO, CO₂, and some other ices from dust grains. These species can linger in the gas phase for a time period much longer than the burst duration, opening a possibility for the indirect detection of the past bursts (Lee 2007; Visser & Bergin 2012; Vorobyov et al. 2013; Jorgensen et al. 2015).

Accretion bursts can also have a strong impact on the structure of the central object depending on the accretion rate and the amount of accretion energy absorbed by the protostar. Preliminary analysis (Baraffe et al. 2009; Baraffe & Chabrier 2010) was based on the assumption of “cold” accretion, assuming that no accretion energy is absorbed by the protostar. In a recent work (Baraffe et al. 2012), we propose a so-called “hybrid” accretion scenario, wherein the amount of accreted energy depends on the accretion rate. This hybrid scenario provides a unified picture linking the observed luminosity spread in young clusters and FU Ori eruptions.

Intense accretion bursts could also account for unexpected lithium and beryllium depletion in some young cluster members (Baraffe & Chabrier 2010; Viallet & Baraffe 2012). Based on arbitrary accretion rates and burst numbers, Baraffe & Chabrier (2010) find that accretion bursts with rates > 10^-4M_⊙ yr^-1 could yield significant Li depletion compared to a non accreting model of same age and mass. A recent observational study from Sergison et al. (2013), based on a careful analysis of star members in the Orion Nebular Cluster and in NGC 2264, confirm the presence of a dispersion in the strength of Li line, as found in previous observational studies. This dispersion indicates a spread in lithium depletion among the members of these clusters, but the level of depletion remains weak. They also find some evidence for weak Li depletion that correlates with luminosity, which could be either explained by an age effect, or by models of past accretion with rates <5×10^-4M_⊙ yr^-1. They however exclude the presence of strongly lithium depleted members in their sample, calling into question the results of Baraffe & Chabrier (2010).

As a step forward to previous models based either on arbitrary or non-consistent accretion rates¹, we present in this paper self-consistent numerical simulations fully coupling numerical hydrodynamics models of collapsing pre-stellar cores (Vorobyov & Basu 2010) and stellar evolution models (Chabrier & Baraffe 1997; Baraffe et al. 1998). We analysed the main impact of consistent accretion history, using hybrid and cold accretion assumptions, on the evolution of young low-mass stars and brown dwarfs and on Li depletion.

2. Model description

We consistently modelled the evolution of the central accreting object and of its circumstellar disk. The model consists of a forming star, described by a one-dimensional (1D) stellar evolution code (see Sect. 2.1), and a circumstellar disk plus infalling envelope described by a two-dimensional (2D) numerical hydrodynamics code (see Sect. 2.2). Both codes are coupled in real time during the simulations (no postprocessing). The input parameter for the stellar evolution code provided by disk modelling is the mass accretion rate onto the star Ṁ. The output of the stellar evolution code is the stellar radius R_∗ and the photospheric luminosity L_∗,ph, which are used by the disk hydrodynamics simulations to calculate the total stellar luminosity and the radiation flux reaching the disk surface. The stellar evolution code is called to update the properties of the protostar every 10 yr of the physical time. For comparison, the global hydrodynamical timestep may be as small as one month.

The coupled evolutions of the disk and of the protostar are followed up to ~1 Myr, when the central object has essentially reached its final mass. Beyond this point we follow the evolution of the latter with the stellar evolution code alone. In some cases, the hydrodynamical simulations predict that the central object is still accreting some material from the disk. In those specific cases, we assume that the mass accretion rate predicted by the last hydrodynamical disk simulation declines linearly to zero during another 1.0 Myr. This is certainly ad-hoc, but follows suggestions from observations of median disk lifetimes ~2−3 Myr with evidence of accretion on the central object still going on (Williams & Cieza 2011). Since at ~1 Myr, the predicted rates from the disk simulations are very small (see Sect. 4), adoption of this arbitrary decay of the rates hardly affects the subsequent evolution of the central object.

2.1. Stellar evolution calculation

The stellar evolution models of accreting protostars are based on the input physics described in Chabrier & Baraffe (1997) and include accretion processes as described in Baraffe et al. (2009 and 2012). The accretion rates onto the protostar are derived from the hydrodynamic calculations described in Sect. 2.2. As in Baraffe et al. (2012), we assume that an amount α of the accretion energy ϵGM_∗Ṁ/R_∗ is absorbed by the protostar, while an amount (1−α) is radiated away and contributes to the accretion luminosity of the star².

In the present calculations, we consider two scenarios as in Baraffe et al. (2012): a cold accretion model with α = 0, meaning that essentially all accretion energy is radiated away, and a hybrid accretion model with α = 0, when accretion rates remain smaller than a critical value Ṁ_cr and α ≠ 0 when Ṁ>Ṁ_cr. More specifically, we adopt $\mathit{\alpha }\mathrm{=}\{\begin{array}{c}\\ \mathrm{0}\mathit{,}& \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{if}\hspace{0.17em}\mathit{Ṁ}\mathit{<}{\mathrm{10}}^{-5}{\mathit{M}}_{\mathrm{\odot }}\mathrm{y}{\mathrm{r}}^{-1}\mathit{,}\\ \mathrm{0.2}\mathit{,}& \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\mathrm{if}\hspace{0.17em}\mathit{Ṁ}\mathrm{\ge }{\mathrm{10}}^{-5}{\mathit{M}}_{\mathrm{\odot }}\mathrm{y}{\mathrm{r}}^{-1}\mathit{.}\end{array}$(1)The justification for a critical value for the accretion rate comes from the analysis of the pressure balance at the stellar surface and is detailed in Appendix B of Baraffe et al. (2012). The choice for the maximum value of α is based on our previous work showing that stellar properties do not change significantly as long as α is greater than 0.1–0.2 (Baraffe et al. 2009, 2012).

Arguments in favour of non zero values for α are provided by the recent radiation hydrodynamics simulations of a protostellar collapse by Vaytet et al. (2013), following the first and second collapse. They find that the accretion shock at the border of the second (Larson) core is subcritical, implying that essentially all the energy provided by the infalling material is absorbed by the core and not radiated away.

For the initial seed mass of the protostar, corresponding to the second Larson core mass, we adopt a value of 1.0M_Jup with an initial radius ~ 1.0 R_⊙. Adopting smaller values for the initial radius would produce even more compact objects (see Sect. 3 and Baraffe et al. 2009, 2012), while choosing a larger radius (> 1 R_⊙) has no significant impact on the evolution. Indeed, the larger the radius, the smaller the thermal timescale (${\mathit{\tau }}_{\mathrm{th}}\mathrm{\propto }\frac{{\mathit{M}}^{\mathrm{2}}}{\mathit{RL}}$) implying that the protostar will rapidly contract and its evolution will be similar to the one starting with 1 R_⊙. We do not vary the initial seed mass, as done in previous works (Hosokawa et al. 2011; Baraffe et al. 2012) on the basis of the calculations of Vaytet et al. (2013) which point toward the universality of the second core mass of about 1 M_Jup.

2.2. Numerical hydrodynamics calculations

The evolution of both the disk and envelope is computed using the numerical hydrodynamics code described in detail in Vorobyov & Basu (2010) and Vorobyov et al. (2013). Here, we briefly recall the main concepts. We start our numerical simulations from the gravitational collapse of a starless cloud core, continue into the embedded phase of star formation, during which a star, disk, and envelope are formed, and terminate our simulations after approximately one Myr of evolution. The protostellar disk, when formed, occupies the inner part of the numerical grid, while the collapsing envelope occupies the rest of the grid. As a result, the disk is not isolated but is exposed to intense mass loading from the envelope. In addition, the mass infall rate onto the disk is not a free parameter of the model but is self-consistently determined by the gas dynamics in the envelope.

To avoid too small time steps, we introduce a sink cell near the coordinate origin with a radius of r_sc = 5 AU and impose a free outflow boundary condition so that the matter is allowed to flow out of the computational domain but is prevented from flowing in. The sink cell is dynamically inactive: it contributes only to the total gravitational potential and ensures a smooth behaviour of the gravity force down to the stellar surface. During the early stages of the core collapse, we monitor the gas surface density in the sink cell and when its value exceeds a critical value for the transition from isothermal to adiabatic evolution, we introduce a central point-mass object representing the forming protostar. From this time on, the stellar evolution code is used to calculate the properties of the protostar. We assume that 90% of the gas that crosses the sink cell lands onto the protostar. A small fraction of this mass (a few per cent) remains in the sink cell to guarantee a smooth transition of the gas surface density across the inner boundary. The other 10% of the accreted gas is assumed to be carried away with protostellar jets. The use of a sink cell implies that the accretion rate inferred at 5 AU is representative of the protostar’s accretion rate. This assumption is a limitation of the models and the reality may be more complex. Additional sources of accretion variability due to thermal instability (Bell & Lin 1994) or magneto-rotational instability (Zhu et al. 2009) may be present at sub-AU scales, changing the properties of the accretion rate.

The main physical processes taken into account when computing the evolution of the disk and envelope include viscous and shock heating, irradiation by the forming star, background irradiation, radiative cooling from the disk surface and self-gravity. Equations and details are provided in Vorobyov & Basu (2010), Vorobyov et al. (2013) and will not be repeated here. The numerical resolution is 256 × 256 grid points and the numerical procedure to solve the hydrodynamics equations is described in detail in Vorobyov & Basu (2010). Irradiation from the central star depends on its luminosity L_∗ which is the sum of the accretion luminosity L_{∗ ,accr} = (1−α)ϵGM_∗Ṁ/R_∗ and the photospheric luminosity L_{∗ ,ph} due to gravitational contraction and eventually deuterium burning in the protostar interior. The stellar mass M_∗ and accretion rate onto the star Ṁ are determined self-consistently during numerical simulations, as described in Vorobyov et al. (2013).

2.3. Initial setup for the hydrodynamical simulations

The initial radial profiles of the gas surface density Σ and angular velocity Ω are based on the derivations from Basu (1997):

$\begin{array}{ccc}& & \\ & & \end{array}$Here, Ω₀ and Σ₀ are the angular velocity and gas surface density at the centre of the core and ${\mathit{r}}_{\mathrm{0}}\mathrm{=}\sqrt{\mathit{A}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathit{\pi G}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{\right)}$ is the radius of the central plateau, where c_s is the initial sound speed in the core. The gas surface density distribution described by Eq. (2) can be obtained (to within a factor of order unity) by integrating the three-dimensional gas density distribution characteristic of Bonnor-Ebert spheres with a positive density-perturbation amplitude A. The value of A is set to 1.2 and the initial gas temperature is set to 10 K.

We have explored a range of initial conditions, characterised by the pre-stellar core mass, the ratio of the rotational to gravitational energy β and the cloud core’s outer radius r_out. We fix r_out/r₀ = 6 in order to generate gravitationally unstable truncated cores of similar form. The adopted values of β lie within the limits inferred by Caselli et al. (2002) for dense molecular cloud cores. To generate the parameters of a specific core with a given value of β, we choose the outer cloud core radius r_out and find r₀ using the adopted ratio above-mentioned. Then, we find the central surface density Σ₀ from the relation ${\mathit{r}}_{\mathrm{0}}\mathrm{=}\sqrt{\mathit{A}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathit{\pi G}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{\right)}$ and determine the resulting cloud core mass M_core from Eq. (2).

3. Parameter survey

The range of initial configuration parameters for all simulations are given in Appendix A for the cold accretion (Table A.1) and hybrid accretion cases (Table A.2). The motivation for the choice of parameters and number of simulations is to produce a population of young stars and brown dwarfs covering the mass range between ~ 0.05 M_⊙ and ~ 1.2 M_⊙ for which observations in various young clusters are available. Because of our assumed density structure for the cores and the truncated condition required for collapse, decreasing the core mass implies increasing the gas surface density Σ₀ and decreasing the radius r₀. Therefore, in order to keep the parameter β within the range of observed values, the value of Ω₀ needs to increase if the core mass decreases. Additionally, in order to get gravitationally unstable disks, we choose values of β close to the upper observational limit, that is ~3%. These choices of parameters and of initial density structure for the cores yield a correlation of increasing Ω₀ with decreasing M_core (see Tables A.1–A.2). We discuss consequences of these choices in Sect. 4.1.

Evolution with time of the mass accretion rates onto the protostar since the beginning of the collapse of the pre-stellar cores for all simulations is given in Figs. A.1 and A.2. To facilitate convergence of the stellar evolution code and avoid too abrupt changes of the accretion rates, actual rates were averaged over the past 50 yr for Ṁ ≤ 10^-5M_⊙ yr^-1 and 25 yr for Ṁ> 10^-5M_⊙ yr^-1. Differences in the time behaviour of Ṁ stem from the different properties of protostellar disks formed from the gravitational collapse of pre-stellar cores. The low-M_core and low-β models produce disks of low mass and size, which are weakly gravitationally unstable and show no sign of fragmentation, while high M_core and β models form disks that are sufficiently massive and extended to develop strong gravitational instability and fragmentation. The forming fragments often migrate onto the star owing to the loss of angular momentum via gravitational interaction with spiral arms or other fragments in the disk, producing strong accretion bursts similar in magnitude to FU-Orionis-type eruptions (e.g. Vorobyov & Basu 2010; Baraffe et al. 2012). In our models, therefore, the stellar accretion and hence the properties of the star are intimately linked to the parental cores and disk properties. Depending on whether accretion is cold (α = 0) or hybrid (α ≠ 0 for Ṁ> 10^-5M_⊙ yr^-1), the impact of accretion on the protostar’s structure is different. As already explained in Baraffe et al. (2009, 2012), cold accretion produces objects more compact and thus less luminous at a given age and mass, compared to both hybrid and non accreting counterparts. Because of the more compact structure of the accreting object in the cold accretion case, the accretion luminosity in this case is larger than in the hybrid case. During the major part of the collapse phase, the accretion luminosity radiated away from the protostar into the disk dominates the photospheric luminosity (Elbakyan et al. 2016). Therefore, despite the lower photospheric luminosity in the cold case, the total luminosity from the protostar radiated into the disc is larger in the cold case during the major duration of the collapse phase. Because of the feedback effect of the protostar luminosity on the disk, consistently taken into account, same initial parameters used in cold and hybrid cases will not produce the same accretion history and the same final object. We however do not find any obvious trend in the property of accretion and burst intensity resulting from the different feedback effects produced by cold and hybrid accretion respectively.

A general test for the assumptions used in the hydrodynamical simulations and for the inferred accretion rates during the embedded phases is provided by the work of Dunham & Vorobyov (2012). It shows that the models provide a reasonable match to the observed protostellar luminosity distribution. We can also compare present results to rates derived from observations (emission lines, accretion shock emission) in young clusters of age ≳1 Myr. This comparison is done in Fig. 1 with model accretion rates predicted at an age of 1 Myr. Observations are based on the recent analysis of the young cluster NGC 2264 (~3−5 Myr) from Venuti et al. (2014). We also show the lower and upper limits of the stellar mass-accretion rate relationship derived for Classical T-Tauri Stars (CTTS) by Fang et al. (2013). Present comparison between models and observations is not a proof of validity of the results, given the large uncertainties both in the hydrodynamical models (see Baraffe et al. 2012) and in the determination of stellar masses, ages and accretion rates from observations. We note that accretion rates for a given mass predicted by the models seem overall larger than the observed rates from Venuti et al. (2014), particularly for low-mass objects. This could be explained by the fact that the model accretion rates are calculated at an age of 1.0 Myr, but the mean age of the NGC 2264 cluster is somewhat older, 3−5 Myr. In addition, as explained in Sect. 3, the low-mass models in our sample have β-values close to the upper observational limit in order for the disks to be gravitationally unstable, which also results in higher (on average) accretion rates. We note that previous comparisons of the model and observed accretion rates for 0.5−3.0 Myr-old objects revealed an overall good agreement (Vorobyov & Basu 2009). It is therefore reassuring to see that essentially all predicted rates are within the lower and upper limits of Fang et al. (2013), and that they follow roughly the observed relations Ṁ ∝ M² in the lower mass range (log M ≲ −0.5) and Ṁ ∝ M in the upper mass range (log M ≳ −0.5) in Fig. 1 (see Fang et al. 2013).

Fig. 1 Relation between accretion rates derived from the hydrodynamical simulations (at an age of 1 Myr) and stellar mass for the cold (blue dots) and hybrid (magenta squares) cases. The magenta crosses are taken from the observations of Venuti et al. (2014). The dashed lines indicate the lower and upper limits observationally derived by Fang et al. (2013).

4. Results

4.1. Analysis of the luminosity spread

In this section we analyse whether the population of synthetic objects resulting from our simulations can naturally produce a luminosity spread in the Herzsprung-Russell diagram (HRD). Figure 2 shows the results for cold and hybrid accretion cases for three different ages, namely 3, 5 and 10 Myr. We exclude models at 1 Myr from the HRD analysis since many of them are still accreting, as illustrated in Fig. 1. Their accretion luminosity is non negligible compared to the photospheric luminosity, and even exceeds it in a few cases. Observationally, most objects in young clusters with inferred ages of ~1 Myr (based on non accreting models, e.g. Baraffe et al. 1998) do not show signatures of strong accretion and their luminosity is essentially dominated by the stellar photosphere thermal emission. This suggests that ages inferred for objects as young as ~1 Myr are very uncertain and strongly model dependent.

In order to illustrate the typical observed luminosity spread, we show in Fig. 2 the data of Bayo et al. (2011, 2012) for the star forming region λ Orionis of typical age ~5 Myr. More thorough comparison of our models with observations in young clusters needs dedicated efforts to obtain observational samples with (i) tested membership; (ii) robust determination of the effective temperature from reliable atmosphere models (presently still lacking) and (iii) case-by-case analysis of very faint objects, which look particularly old, to eliminate objects with edge-on disks (as done in e.g. Sergison et al. 2013).

Fig. 2 Position of cold and hybrid models in a HRD at given ages. The black dashed lines provide isochrones for stellar ages of 1 Myr, 3 Myr, 5 Myr, 10 Myr and 50 Myr (from top to bottom) derived from non-accreting stellar evolution models of Baraffe et al. (1998). The cross symbols are observations from Bayo et al. (2011, 2012) in λ Orionis (~5 Myr).

Two striking features are displayed in Fig. 2. The first is that instead of showing a spread, the models tend to show a bimodal distribution: hybrid models at a given age tend to cluster along the non-accreting isochrone of same age. The cold accretion models at a given age are systematically fainter and look significantly older (5 to 10 times older) compared to the non-accreting isochrone of same age. The second feature is the lack of luminosity dispersion for log T_eff ≲ 3.5, that is T_eff≲ 3200 K. This corresponds to masses ≲0.1 M_⊙. Cold and hybrid accretion produce essentially the same objects in the very low mass regime because of the low accretion variability predicted by the simulations and the imposed value of Ṁ_cr which strongly limits the number of bursts during which the accreting object can absorb some amount of the accretion energy (see Fig. A.2).

One can fill the gaps by varying one parameter in the hybrid scenario, namely the critical accretion rate Ṁ_cr above which accretion switches from cold to hot. We recall that fiducial values in our hybrid scenario are Ṁ_cr = 10^-5 and α = 0.20. To test the sensitivity of the results to this parameter, we have rerun two fully coupled sequences with parameters of models 10 and 28, respectively, but with Ṁ_cr = 5 × 10^-5 and one sequence with parameters of model 6 but with Ṁ_cr = 3 × 10^-6 (see Table A.2 for the corresponding parameters). Increasing the value of Ṁ_cr limits the heating of the protostar due to the absorption of accretion energy during bursts, and thus produces an object with a structure intermediate between the ones produced by cold and fiducial hybrid scenarios, respectively. This helps filling the gap in the HRD for objects with log T_eff ≳ 3.5 as illustrated in Fig. 3. We note that decreasing the value of α would have the same effect. On the opposite, decreasing the value of Ṁ_cr allows for more absorption of accretion energy, and thus more heating of the structure of the accreting object at the low end of the mass distribution, limiting the effect of mass accretion and producing less compact objects. This helps producing a luminosity spread for the lowest mass objects with log T_eff ≲ 3.5 as illustrated in Fig. 3. Justifications for varying these two parameters α and Ṁ_cr are discussed in Sect. 5.

As mentioned in Sect. 3, our choice of parameters for the initial core properties favours fast rotating low mass cores. More slowly rotating cores, however, certainly exist. They will form low mass objects with small stable disks, with steady instead of burst-like accretion. If, under such conditions, the accretion timescale, ${\mathit{\tau }}_{\mathrm{acc}}\mathrm{=}\frac{\mathit{M}}{\mathit{Ṁ}}$, remains smaller than the thermal timescale, τ_th, of the protostar, these objects will be representative of the cold accretion case and explain the faintest observed young cluster members (see Fig. 2 and Baraffe et al. 2009, 2012). In the opposite case, τ_acc > τ_th, the objects will have a larger radius, thus a larger luminosity. Therefore, this population of small slowly rotating cores, not explored in the present study, will help resolving the problem of reduced luminosity spread in the low-mass domain.

Fig. 3 Effect of varying the value of Ṁcr in the hybrid accretion scenario. The curves are the same as Fig. 2.

Fig. 4 Surface lithium abundance (in units of the initial mass fraction Li = 10-9) as a function of effective temperature in accreting models under the cold (left panels) and hybrid (right panels) accretion scenarios. The dashed curves are the predictions from non-accreting stellar evolution models of Baraffe et al. (1998). The latter is not visible in the 3 Myr panels since non-accreting models do not predict any depletion at this age and thus Li/Li0 = 1 for all models.

Fig. 5 Same as Fig. 4 but as a function of mass.

Fig. 6 Evolution of various quantities as a function of time (in yr) for model 29 (dash-dotted magenta curve) with final mass 0.735 M⊙ and model 30 (solid blue curve) with final mass 0.752 M⊙ under the assumption of cold accretion. Li/Li0 is the surface abundance of lithium in units of the initial mass fraction. Mrad is the size of the radiative core in M⊙. Tc is the central temperature and Ṁ (in M⊙ yr-1) is the accretion rate predicted by the hydrodynamical simulations. The black long-dashed curves are the predictions from a non-accreting stellar evolution models with mass 0.75 M⊙ from Baraffe et al. (1998).

4.2. Lithium depletion

Present consistent calculations enable to re-analyse on more solid grounds the results of Baraffe & Chabrier (2010) about lithium depletion. The possibility of having strong lithium depletion due to early accretion needs to be confirmed based on more consistent model assumptions since this prediction was called into question by the observations of Sergison et al. (2013). Results are displayed in Figs. 4 and 5 which show the surface abundance of lithium as a function of effective temperatures and stellar mass, respectively, for three different ages. We summarise the main results below.

• (a) Both cold and hybrid accretion scenarios confirm thepossibility of accelerating lithium depletion at ages of afew Myr compared to non accreting models. Asexplained in Baraffe & Chabrier (2010), thisstems from the more compact and hotter structure produced bymass accretion. Hybrid accretion, as expected, tends to moderatethis effect, producing less anomalously Li-depleted objectscompared to non-accreting models.

• (b) Among all calculated models, we find only one extremely depleted case compared to the non accreting counterpart (Li/Li₀< 1% at t~ 1 Myr). This is model 29 in the cold case (final mass 0.735 M_⊙)³. Interestingly enough, this sequence is not characterised by strong accretion bursts, as shown in the lower panel of Fig. 6 (see point (c) below).

• (c) Baraffe & Chabrier (2010) obtained extremely high level of depletion when the accreting object undergoes from the beginning of its evolution consecutive violent bursts exceeding ~10^-4M_⊙ yr^-1. The hydrodynamical disk evolution models predict the existence of strong bursts in excess of 10^-4M_⊙ yr^-1. However, models based on consistent accretion rates show no clear correlation between the existence of such strong bursts and anomalous level of Li depletion. For example, model 25 (final mass 0.53 M_⊙) with cold accretion displays several extreme bursts up to ~10^-3M_⊙ yr^-1 (see Fig. A.1) but shows no abnormally high Li depletion at an age of a few Myr (Li/Li₀ = 0.993 at 3 Myr). On the other hand, model 33 (final mass 0.926 M_⊙) with cold accretion shows similar intense bursts and has anomalously depleted its Li by more than a factor three at 3 Myr (Li/Li₀ = 0.286 at 3 Myr). In the case of hybrid accretion, as mentioned in (a) above, models with the strongest bursts show no anomalous Li depletion.

• (d) As a new result, Figs. 4 and 5 show accreting models with higher surface Li abundance than non-accreting counterparts at same age. Those are models in the mass range ~0.5–1.1 M_⊙ (see Fig. 5) which develop a radiative core during their evolution. This unexpected result is due to a subtle competition between different effects: (i) the more compact structure and the consequent increase of central temperature due to mass accretion; (ii) the growth of the radiative core which strongly depends on the central temperature and density (see e.g. Baraffe & Chabrier 2010) and (iii) the accretion of fresh lithium from the accretion disk.

• (e) The results show that all models with higher Li depletion compared to their non accreting counter-part appear to be fainter in a HRD and thus look older than the non-accreting models. But we stress that there is no systematic correlation between the position in the HRD and the level of lithium depletion, meaning that all objects that look “older” are not necessarily anomalously depleted in Li. As an example, the cold accretion models 22 (final mass 0.448) and 23 (final mass 0.470) have the same position in the HRD, lying at 5 Myr on the 50 Myr non-accreting isochrone. But model 22 has a Li surface abundance depleted by a factor ~2 whereas model 23 shows essentially no Li depletion at 5 Myr.

This ensemble of results shows that there is no simple rule or trend linking the intensity of an accretion burst, the level of Li depletion and the position in the HRD. The evolution of Li is very sensitive to the strength of an accretion burst and whether the burst happens when a radiative core is already developed. But it also depends on how rapidly the convective envelope shrinks and on the relative amount of fresh Li that will still be accreted compared to the abundance of Li in the shrinking convective envelope. Such high sensitivity of results is illustrated in Fig. 6: at time t~ 0.1 Myr model 30 undergoes several strong accretion bursts which rapidly heat up the central stellar region and cause a rapid growth of the radiative core before a significant amount of Li could be depleted. The additional accretion of fresh Li in a convective envelope which is smaller than that of model 29 or of the non accreting counterpart, yields a young object with a Li abundance in the convective envelope that is essentially the same as the initial Li abundance in the collapsing pre-stellar core.

5. Discussion and conclusion

The self-consistent evolutionary models of accreting objects presented in this work⁴ confirm the prediction of a spread in luminosity in HRDs for a coeval population of young low-mass stars and brown dwarfs and the possibility to accelerate lithium depletion. Admittedly, the spread in the HRD predicted by present models is not fully satisfactory. It resembles a bimodal distribution rather than a true spread in luminosity. Cold accretion models are systematically fainter than the non-accreting isochrone of same age while hybrid accretion models lie close to it. Limiting the values of free parameters in the present self-consistent calculations with α fixed either to zero (cold) or to 0.2 (hybrid), and a fixed value for Ṁ_cr is too simplistic and only provides limiting cases. The fact that, within this restricted choice of parameters, the models produce a limited luminosity spread could be solved by using different values for α and Ṁ_cr.

Determination of the value of α characterising the amount of energy transferred to a central object by material accreted from a disk is a formidable radiative hydrodynamical problem that is far from being solved. In the context of young low mass objects, whether accretion is cold or hot is an open debate. Currently, neither models nor observations can solve this question (see e.g. Hartmann et al. 2011; Baraffe et al. 2012).

Assuming that accreted material does reach the surface of the central object with some amount of internal energy, justifications for varying Ṁ_cr, which delineates the regimes between cold and hot accretion, can be found in the recent multi-dimensional simulations of accretion onto young stars by Geroux et al. (2016). First, these multi-D simulations provide a test for the assumptions used in 1D evolutionary calculations of accreting objects as presented in this work. They are reassuring regarding the validity of qualitative results based on 1D accreting models. One interesting finding is that the treatment of hot accretion (for α≳ 0.1) is likely more realistic if based on an outer accretion boundary condition at the stellar surface rather than assuming redistribution of accretion energy deep in the interior, as presently assumed in our 1D evolutionary calculations. Geroux et al. (2016) show that for a given accretion rate and a given value of α, the later assumption may overestimate the effect on the structure, in terms of radius expansion of the object. But the important conclusion is that the effects described by the multi-D simulations for a given α can be mimicked by a treatment assuming redistribution of energy, provided a smaller value of α is used.

The multi-D simulations also show that there is a change in behaviour at a given value of α in the redistribution of accreted mass and energy in the stellar interior and in the effect of accretion on the stellar structure. The threshold value of α characterising this change in behaviour depends on the mass accretion rate and on the difference between the entropy of the accreted material and the bulk entropy in the convective envelope of the accreting object. It thus highly depends on the structure of the accreting object, and thus on its mass. This suggests that the effect of hot accretion, which essentially expands the object and compensates for the effect of mass accretion (which makes the object more compact and fainter) cannot be mimicked by one single value of α and a fixed value of Ṁ_cr, but more likely by values that vary depending on the accretion rate and the mass (or the interior bulk entropy) of the accreting object. Work is in progress to explore the implications of the multi-D simulation results. The goal is to examine whether the use of an accretion boundary in the 1D evolutionary calculations can more naturally produce a luminosity spread than using the assumption of redistribution of accreted energy in the stellar interior and varying both the values of α and Ṁ_cr.

The analysis of Li depletion suggests that extreme Li depletion due to early accretion should be a very rare event. Among ~60 models calculated under various assumptions and parameters, we only find one model which displays almost complete Li depletion. This is not in contradiction with the results of Sergison et al. (2013) who find no strong Li depleted objects among ~168 targets in NGC 2264 and the Orion Nebular Cluster. Interestingly enough, Bouvier et al. (2016) recently re-analysed Li abundance in NGC 2264 members with a sample of ~200 stars. They report one object with low lithium equivalent width but with radial velocity consistent with the cluster’s average and with rotational velocity vsini consistent with a young age. This source was rejected from the Bouvier et al. (2016) analysis but we urge its follow-up with further membership analysis, since confirmation of its membership could provide a genuine signature of our accretion scenario. Our results also show that models with high Li depletion compared to non-accreting models are also fainter in a HRD. But we stress that all objects that look older in a young cluster should not necessarily show a lower level of Li depletion.

The interesting result about Li depletion in accreting models is the prediction of Li-excess compared to surface Li abundance in the non-accreting counterparts. This excess is due to a subtle competition between the effect of burst accretion and its impact on the central temperature, the growth of the radiative core and the accretion of fresh Li. Only consistent models could reveal such a subtle combination of effects. In their analysis of NGC 2264, Bouvier et al. (2016) suggest a connection between lithium abundance and rotation, with Li-excess in fast rotators compared to slow ones. This discovery is at odd with current understanding of rotation in stars with rotational mixing enhancing Li depletion. Our results suggest a possible explanation for this puzzle. Strong accretion bursts that occur when the central object develops a radiative core could both explain the higher spin of the star, as a result of angular momentum accretion, and its slower Li depletion due to rapid shrinkage of the convective envelope and accretion of fresh Li. This idea is attractive and deserves further exploration by both observations and models.

To finish, both hydrodynamical and evolutionary models presented in this work contain many uncertainties, part of them already discussed in Baraffe et al. (2012). Recent multi-D simulations of realistic stellar structures as presented in Geroux et al. (2016) provide some preliminary tests of major assumptions used in 1D stellar evolution models of accreting objects and give some more confidence in their qualitative predictions. Much more work needs to be done before we can believe in quantitative results and build realistic synthetic populations of young objects accounting for their early accretion history to compare with observations. We believe this work provides one step forward in this direction. But guidance from observations are key to advance on the front of model developments. We suggest within

this context more observational efforts devoted to characterising the properties of young objects undergoing accretion bursts such as FUors and EXors. Lastly, the standard practice of rejecting outliers with abnormal Li depletion in young clusters may be a lost opportunity to find a genuine confirmation of our accretion scenario. We strongly advise for a follow-up studies of such outliers.

Acknowledgments

We thanks Amelia Bayo, Min Fang and Laura Venuti for providing their data. This project was partly supported by the European Research Council through grants ERC-AdG No. 320478-TOFU and No. 247060-PEPS, by the Russian Ministry of Education and Science Grant 3.961.2014/K and by the Austrian Science Fund (FWF) under research grant I2549-N27. V.G.E. acknowledges Southern Federal University Development Program for financial support. The simulations were performed on the Vienna Scientific Cluster (VSC-2), on the Shared Hierarchical Academic Research Computing Network (SHARCNET), on the Atlantic Computational Excellence Network (ACEnet) and on the University of Exeter Supercomputer, a DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS and the University of Exeter.

References

Appendix A: Initial model data

We provide the details of the initial parameters of pre-stellar cores in cold accretion (Table A.1) and in the hybrid accretion simulations (Table A.2). In both tables, the second column is the initial core mass, third column is the ratio of rotational to gravitational energy, fourth column is the central angular velocity, fifth column is the radius of the central near-constant-density plateau, and sixth column is the central surface density. In addition, the seventh column provides the final stellar masses in each model. The evolution with time of the mass accretion rates onto the protostar since the beginning of the collapse of the pre-stellar cores for all simulations are shown in Figs. A.1 and A.2.

Fig. A.1 Mass accretion rates vs. time for the 36 models calculated for cold accretion. The horizontal dotted lines mark the critical value for the transition from cold to hot accretion.

Fig. A.2 Mass accretion rates vs. time for the 31 models calculated with hybrid accretion. The horizontal dashed lines mark the critical value for the transition from cold to hot accretion.

All Tables

All Figures

	Fig. 1 Relation between accretion rates derived from the hydrodynamical simulations (at an age of 1 Myr) and stellar mass for the cold (blue dots) and hybrid (magenta squares) cases. The magenta crosses are taken from the observations of Venuti et al. (2014). The dashed lines indicate the lower and upper limits observationally derived by Fang et al. (2013).
In the text

	Fig. 2 Position of cold and hybrid models in a HRD at given ages. The black dashed lines provide isochrones for stellar ages of 1 Myr, 3 Myr, 5 Myr, 10 Myr and 50 Myr (from top to bottom) derived from non-accreting stellar evolution models of Baraffe et al. (1998). The cross symbols are observations from Bayo et al. (2011, 2012) in λ Orionis (~5 Myr).
In the text

	Fig. 3 Effect of varying the value of Ṁcr in the hybrid accretion scenario. The curves are the same as Fig. 2.
In the text

Fig. 4 Surface lithium abundance (in units of the initial mass fraction Li = 10-9) as a function of effective temperature in accreting models under the cold (left panels) and hybrid (right panels) accretion scenarios. The dashed curves are the predictions from non-accreting stellar evolution models of Baraffe et al. (1998). The latter is not visible in the 3 Myr panels since non-accreting models do not predict any depletion at this age and thus Li/Li0 = 1 for all models.

In the text

	Fig. 5 Same as Fig. 4 but as a function of mass.
In the text

Fig. 6 Evolution of various quantities as a function of time (in yr) for model 29 (dash-dotted magenta curve) with final mass 0.735 M⊙ and model 30 (solid blue curve) with final mass 0.752 M⊙ under the assumption of cold accretion. Li/Li0 is the surface abundance of lithium in units of the initial mass fraction. Mrad is the size of the radiative core in M⊙. Tc is the central temperature and Ṁ (in M⊙ yr-1) is the accretion rate predicted by the hydrodynamical simulations. The black long-dashed curves are the predictions from a non-accreting stellar evolution models with mass 0.75 M⊙ from Baraffe et al. (1998).

In the text

	Fig. A.1 Mass accretion rates vs. time for the 36 models calculated for cold accretion. The horizontal dotted lines mark the critical value for the transition from cold to hot accretion.
In the text

	Fig. A.2 Mass accretion rates vs. time for the 31 models calculated with hybrid accretion. The horizontal dashed lines mark the critical value for the transition from cold to hot accretion.
In the text