© ESO, 2017

1. Introduction

The evolution of the pre-main-sequence (PMS) stars has been studied for decades by various authors (e.g., Baraffe et al. 1998; Palla & Stahler 2000). Regardless of some model uncertainties, stellar evolution calculations show consensus for the standard evolution path, which extends through the Hayashi and Henyey tracks on the Hertzsprung-Russell (HR) diagram (Henyey et al. 1955; Hayashi 1961). Effects of mass accretion have also been incorporated in stellar evolution calculations (e.g., Stahler et al. 1980; Palla & Stahler 1991; Hosokawa & Omukai 2009). Palla & Stahler (1990) proposed the concept of the protostellar birth line, where the PMS stars first become optically visible as the surrounding envelope disperses and the mass accretion ceases.

Recent studies further added updates considering the mass accretion in more realistic situations. First of all, while most of the previous studies assumed simple constant accretion histories for simplicity, numerical simulations are revealing a more sporadic nature of the mass accretion. Since parental molecular cores have finite amounts of angular momentum, a protostar accretes the gas through a circumstellar disk that forms thanks to the near-conservation of the angular momentum in the collapsing core. In the early embedded phase of disk evolution, where most of the final stellar mass is accumulated, the angular momentum and mass transport is mostly driven by the gravitational torque (Vorobyov & Basu 2009b). In this phase, the disks are often prone to gravitational instability, which results in very time-dependent mass accretion histories, for example, the episodic accretion histories that repeat short accretion burst events and relatively longer quiescent phases (Vorobyov & Basu 2010, 2015; Machida et al. 2011; Tomida et al. 2017). Other mechanisms, such as the magnetorotational and thermal instabilities, planet-disk interactions, and close stellar encounters, can also produce episodic accretion bursts in both the embedded and T Tauri phases of disk evolution (see Audard et al. 2014, for a review). There are also various observational signatures of such episodic accretion reported for low-mass protostars (e.g., Dunham et al. 2010; Liu et al. 2016).

Effects of such variable mass accretion have recently been included in stellar evolution calculations that focus on the early (Vorobyov et al. 2017; Hosokawa et al. 2016) and late evolutionary phases (e.g., Baraffe et al. 2009, 2012, 2016; Hosokawa et al. 2011). These studies reported that PMS tracks can largely differ from the non-accreting tracks that have been used to estimate the stellar age spreads in young clusters (e.g., Soderblom et al. 2014). In some cases, a PMS star first appears quite near the main sequence on the HR diagram immediately after the mass accretion ceases, while the non-accreting isochrones assume that the PMS star starts the evolution in an upper part on the HR diagram. As a result, the age of the star as inferred from the non-accreting tracks can be overestimated by about a factor of 2.

Although the above mismatch was found in studies with variable mass accretion, another factor is normally incorporated with a free parameter in the model: the thermal efficiency of the accretion, that is, the specific entropy of the accreted gas. It has been suggested that this “warmth” of the accretion controls the PMS evolution instead of the time-variability of the accretion rates (e.g., Hartmann et al. 1997, 2016; Hosokawa et al. 2011; Kunitimo et al. 2017). However, these effects are often mixed because the accretion thermal efficiency can be modeled as a function of the accretion rates (Baraffe et al. 2012). It is still quantitatively uncertain to which extent the non-accreting isochrones are reliable after all, and whether the accretion variability or thermal efficiency of the accretion plays the main role in the PMS evolution.

In this work, we revisit the PMS evolution of (sub-)solar mass stars and brown dwarfs by taking variable mass accretion into consideration and making different assumptions regarding the thermal efficiency parameter α, the fraction of the accretion energy taken into the stellar interior. We use variable mass accretion histories on a protostar (Ṁ) obtained in numerical simulations following the formation and evolution of protostars and protostellar disks described in detail in Vorobyov & Basu (2010). Then we compute the evolution of the protostars using a modified version of the stellar evolution code STELLAR (Yorke & Bodenheimer 2008; Hosokawa et al. 2013; Sakurai et al. 2015). This approach is similar to that employed in Baraffe et al. (2012), but using the STELLAR code rather than the Lyon code (Baraffe et al. 2009). We compare the PMS evolutionary tracks with the above variable accretion to the non-accreting isochrones calculated with the same code STELLAR, and monitor how large the offsets are that the accreting tracks show with the isochrones and isomasses as functions of the time. Moreover, we also quantitatively analyze how this evolution of the offsets varies with different modeling of the thermal efficiency α. Our results show that with α> 10^-3 (we call these hybrid and hot accretion cases below), offsets with the isochrones and isomasses only appear for early ages of <a few Myr, where even the non-accreting isochrones suffers from relatively large modeling uncertainties. At later ages, the notable offsets only remain when α is fixed at 10^-3 (we call this cold accretion below), clearly showing that the value of α is a key for altering the PMS tracks.

The rest of the paper is organized as follows. We first explain our adopted methods on the numerical simulations of the disk accretion and on the stellar evolution calculations in Sects. 2 and 3, respectively. The results are described in Sects. 4–6, where the cases are divided in terms of the different thermal efficiencies, that is, the hybrid, hot, and cold accretion. In Sect. 7 we propose possible updates of the isochrones considering the effects of mass accretion. Finally, Sects. 8 and 9 are devoted to the discussions and conclusions.

2. Model description

We computed the evolution of low-mass stars and upper-mass brown dwarfs¹ using the stellar evolution code STELLAR that was originally developed by Yorke & Bodenheimer (2008), with recent modifications as described in Sect. 3. We started our computations from a protostellar seed, continued through the main accretion phase where the growing star accumulates most of its final mass, and completed our computations when the star approached the main sequence. We used the protostellar accretion rates obtained from the numerical hydrodynamics simulations of prestellar core collapse described in Sect. 2.1. We note that this approach is not fully self-consistent. The stellar evolution calculations are not fully coupled with the disk evolution, but employ the precomputed mass accretion rates. The full coupling of the disk evolution modeling with stellar evolution calculations was recently reported in Baraffe et al. (2016) using the Lyon code. Such real-time coupling has also been performed with the STELLAR code, mostly to study the high-mass star formation (e.g., Kuiper & Yorke 2013; Hosokawa et al. 2016). We also plan to apply the STELLAR code in future studies of the low-mass star formation.

2.1. Numerical hydrodynamics code

The evolution of protostellar disks was computed using the numerical hydrodynamics code described in detail in Vorobyov & Basu (2010). Here, we briefly review the main concepts and equations. We start our numerical simulations from the gravitational collapse of a gravitationally contracting cloud core, continue into the embedded phase of star formation, during which a star, disk, and envelope are formed, and terminate our simulations after 1.0–2.0 Myr of evolution, depending on the model and available computational resources². 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 the embedded phase of star formation. 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 mass accretion rate is calculated as the mass passing through the sink cell per time step of numerical integration. 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 main physical processes taken into account when computing the evolution of the disk and envelope include viscous and shock heating, irradiation by the forming star, background irradiation, radiative cooling from the disk surface, and self-gravity. The code is written in the thin-disk limit, complemented by a calculation of the disk vertical scale height using an assumption of local hydrostatic equilibrium. The resulting model has a flared structure, which guaranties that both the disk and envelope receive a fraction of the irradiation energy from the central protostar. The corresponding equations of mass, momentum, and energy transport are where subscripts p and p′ refer to the planar components (r,φ) in polar coordinates, Σ is the mass surface density, e is the internal energy per surface area, P is the vertically integrated gas pressure calculated via the ideal equation of state as P = (γ−1)e with γ = 7/5, is the velocity in the disk plane, and is the gradient along the planar coordinates of the disk. The gravitational acceleration in the disk plane, , takes self-gravity of the disk into account, which is found by solving for the Poisson integral (see details in Vorobyov & Basu 2010), and the gravity of the central protostar when formed. Turbulent viscosity due to sources other than gravity is taken into account via the viscous stress tensor Π, the expression for which is provided in Vorobyov & Basu (2010). We parameterize the magnitude of kinematic viscosity ν using the alpha prescription with a spatially and temporally uniform α_visk = 5 × 10^-3.

The radiative cooling Λ in Eq. (3) is determined using the diffusion approximation of the vertical radiation transport in a one-zone model of the vertical disk structure, (4)where σ is the Stefan-Boltzmann constant, T_mp = Pμ/RΣ is the midplane temperature of gas, μ = 2.33 is the mean molecular weight (90% of atomic H, 10% of atomic He), R is the universal gas constant, and F_c = 2 + 20tan^-1(τ)/(3π) is a function that secures a correct transition between the optically thick and optically thin regimes. We use the frequency-integrated opacities of Bell & Lin (1994). The heating function is expressed as (5)where T_irr is the irradiation temperature at the disk surface determined by the stellar and background blackbody irradiation as (6)where T_bg is the uniform background temperature (in our model set to the initial temperature of the natal cloud core) and F_irr(r) is the radiation flux (energy per unit time per unit surface area) absorbed by the disk surface at radial distance r from the central star. The latter quantity is calculated as (7)where L_∗ is the total (accretion plus photospheric) luminosity of the central object (star or brown dwarf), and γ_irr is the incidence angle of radiation arriving at the disk surface at a radial distance r. The incidence angle is calculated using the disk surface curvature inferred from the radial profile of the disk vertical scale height (see Vorobyov & Basu 2010, for more details). The photospheric luminosity L_ph is calculated from the D’Antona & Mazzitelly stellar evolution tracks (D’Antona & Mazzitelli 1994), while the accretion luminosity is calculated as follows: (8)where M_∗ is the mass of the central object and R_∗ the radius of the central object (also provided by the D’Antona & Mazzitelly tracks). We note that D’Antona & Mazitelli’s tracks do not cover the very early phases of stellar evolution. Therefore, we have used a power-law expression to extrapolate to times earlier than those included in the pre-main-sequence tracks, L_ph = L_ph,0(t/t0)⁴, where t₀ is the earliest time in the tracks and L_ph,0 is the pre-main-sequence luminosity at this time. We assume ϵ = 1/2, which is characteristic of accretion from a thin disk. Because we use the non-accretion isochrones when calculating the mass accretion histories, the thermal efficiency of accretion α cannot be taken into account. However, in the stellar evolution calculations using the precomputed accretion histories, we do take the thermal efficiency of accretion into account (see Sect. 3). Equations (1)–(3) are solved on the polar grid (r,φ) with a numerical resolution of 512 × 512 grid points. The numerical procedure is described in detail in Vorobyov & Basu (2010).

2.2. Initial conditions for the numerical hydrodynamics code

The initial radial profiles of the gas surface density Σ and angular velocity Ω of collapsing prestellar cores have the following form: Here, Σ₀ and Ω₀ are the gas surface density and angular velocity at the center of the core and is the radius of the central plateau, where c_s is the initial sound speed in the core. These radial profiles are typical of prestellar cores that are formed as a result of the slow expulsion of the magnetic field due to ambipolar diffusion, with the angular momentum remaining constant during axially symmetric core compression (Basu 1997). The value of the positive density perturbation A is set to 1.2, and the initial gas temperature in collapsing cores is set to 10 K.

Each model is characterized by a distinct ratio r_out/r₀ = 6 in order to generate gravitationally unstable truncated cores of similar form, where r_out is the cloud core’s outer radius. The actual procedure for generating the parameters of a specific core consists of two steps. First, we choose the outer cloud core radius r_out and find r₀ using the adopted ratio between these two quantities. Then, we find the central surface density Σ₀ from the relation and determine the resulting cloud core mass M_cl from Eq. (9). When the gas surface density profile and the initial mass are fixed, we set the angular velocity profile by choosing the value of Ω₀ so that the ratio of rotational to gravitational energy β falls within the limit inferred by Caselli et al. (2002) for dense molecular cloud cores. The parameters of our models are presented in Table 1. The second column is the initial core mass, the third column is the ratio of rotational to gravitational energy, the fourth column is the radius of the central near-constant-density plateau, and the fifth column is the final stellar mass. The models are ordered in the sequence of increasing final stellar masses M_{∗ , fin}.

2.3. Accretion rate histories

The available computational resources allow us to calculate the disk evolution, and hence, the protostellar accretion history only for about 1.0–2.0 Myr. We assume that in the subsequent evolution the mass accretion rate declines linearly to zero during another 1.0 Myr. Effectively, this implies a disk dispersal time of 1 Myr and a total disk age of about 2.0–3.0 Myr in our models. These values are in general agreement with the disk ages inferred from observations (Williams & Cieza 2011), although some objects demonstrate longer disk lifetimes.

We calculate protostellar accretion histories for the 35 prestellar cores listed in Table 1. Figure 1 presents Ṁ vs. time elapsed since the beginning of numerical simulations. The horizontal dashed lines mark a critical value of the mass accretion rate Ṁ_cr = 10^-5M_⊙ yr^-1 above which accretion changes from cold to hot in the hybrid accretion scenario (see Sect. 3.2). The protostellar accretion rate is calculated as the mass passing through the sink cell per unit time, Ṁ = −2πr_scΣv_r, where v_r is the radial component of the velocity at the inner computational boundary.

Fig. 1 Mass accretion rates vs. time for 35 models. The horizontal dashed lines mark the critical value for the transition from cold to hot accretion scenario. The models are arranged so that the final stellar mass increases from left to right and from top to bottom. The vertical dashed, dotted, and dash-dotted lines mark time instances when 90%, 95%, and 98% of the total stellar mass has been accumulated.

Evidently, the accretion rates in the low-M_core and low-β models (e.g., models 8, 11, 13, 18, and 21) are characterized by an order-of-magnitude flickering that gradually diminishes with time. On the other hand, models with high M_core and β (e.g., models 14, 16, 22, and 23) demonstrate large-amplitude variations in Ṁ and strong accretion bursts exceeding in magnitude 10^-5M_⊙ yr^-1. This difference in the time behavior of Ṁ stems from the different properties of protostellar disks formed from the gravitational collapse of prestellar cores (Vorobyov 2010). 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 (Vorobyov & Basu 2006, 2010, 2015).

Fig. 2 Stellar evolution sequences on the total luminosity L∗ – effective temperature Teff diagram for the hybrid accretion models. The dots present the model tracks and the color of the dots varies according to the stellar age shown in the vertical bar (in log yr). The zero-point age for each model is marked by the red diamonds. The green symbols mark the reference ages, as indicated in the bottom left corner, that has elapsed since the zero-point age of each object. The black solid lines provide the isochrones for stellar ages of 1 Myr, 10 Myr, and 25 Myr (from top to bottom) derived from the non-accreting stellar evolution models of Yorke & Bodenheimer (2008).

3. Stellar evolution code

We used the stellar evolution code STELLAR originally developed by Yorke & Bodenheimer (2008). The detailed description of the code can be found in Sakurai et al. (2015). We here briefly review the main features of the code.

The basic equations to be solved are as follows: where m is the mass contained within a spherical layer with radius r, P is the total (gas plus radiation) pressure, l is the local luminosity, E_nut is the specific energy production rate by nuclear reactions, c_P is the isobaric specific heat, T is the temperature, δ ≡ −(∂lnρ/∂lnT)_P, and ∇ ≡ ∂lnT/∂lnP is the temperature gradient calculated using the mixing-length theory for convective layers. Nuclear reactions are computed up to the helium burning (3α and {CNO}+He).

The model stars are divided into two parts: the atmosphere, which provides surface boundary conditions, and the interior, which contains most of the stellar material. We consider the atmosphere to be spherically symmetric and gray, applying the equations of hydrostatic equilibrium and radiative or convective energy transport, which are integrated inward by a Runge-Kutta method. We integrate the atmospheric structure down to a fitting point, where the boundary conditions for the interior are provided. We use the standard Henyey method to cause the interior structure to converge to satisfy the boundary conditions. As the accreting gas settles onto the stellar surface, the added material gradually sinks inward through the atmosphere and is eventually incorporated into the interior by a rezoning procedure.

3.1. Initial conditions for the stellar evolution code

The stellar evolution calculations started from a fully convective polytropic stellar seed of 5 Jupiter masses and 3 Jovian radii with a polytropic index n = 1.5. Before commencing actual calculations that take the mass growth via accretion into account, the polytropic seed was allowed to relax to a fully converged stellar model. The composition of the gas was assumed to be X = 0.71 and Y = 0.273. The deuterium abundance was set to 40 ppm. The initial effective temperature of the seed was 1200 K.

3.2. Thermal efficiency of accretion

After setting the properties of the initial protostellar seed, we commenced the stellar evolution calculations using the accretion rate histories calculated in Sect. 2.3. During these calculations, we assumed that a fraction α of the accretion energy ϵGM_∗Ṁ/R_∗ is absorbed by the protostar, while a fraction 1−α is radiated away and contributes to the accretion luminosity of the star. Here, M_∗ and R_∗ are the mass and radius of the central star. In this paper, we consider three scenarios for the thermal efficiency of accretion: (i) cold accretion with a constant α = 10^-3, meaning that practically all accretion energy is radiated away and little is absorbed by the star, (ii) hot accretion with a constant α = 0.1, and (iii) a hybrid scheme defined as follows: (15)This functional form of α has the property that accretion remains cold at small Ṁ and gradually changes to hot accretion above a certain critical value of Ṁ_cr, for which we chose 10^-5M_⊙ yr^-1 based on modeling of FU Orionis-type eruptions of Kley & Lin (1996) and Hartmann et al. (2011) and analytical calculations of Baraffe et al. (2012), which showed that such a transition reflects a change in the accretion geometry with increasing accretion rate and/or in the magnetospheric interaction between the star and the disk. We note that the actual value of α in the hybrid accretion scheme changes smoothly from 10^-3 to 0.1 over a time period of ~1.0 Myr to preserve numerical stability of the stellar evolution code.

4. Hybrid accretion

In this section, we present our model results for the hybrid accretion scenario, in which the value of α depends on the mass accretion rate. The computed stellar evolution sequences of the total (accretion plus photospheric) luminosity L_∗ vs. effective temperature T_eff for 31 models³ are shown with the colored dots in Fig. 2. The color of the dots varies according to the stellar age that is indicated on the vertical bar (in log yr). The zero-point of the stellar age is defined as the time instance when the growing star accumulates 95% of its final mass. At this time instance, the envelope is essentially dissipated and the disk-to-star mass ratio is 5%. This value is in agreement with the inferred upper limit on the disk-to-star mass ratios of young T Tauri stars in nearby star-forming regions (e.g., Palla & Stahler 2000), especially considering that the disk masses may in fact be somewhat underestimated when using dust continuum measurements (Dunham et al. 2014; Tsukamoto et al. 2016). In other words, the adopted 95% criterion for the zero age corresponds to young T Tauri stars that have just left the protostellar phase of stellar evolution. We also emphasize the qualitative change in character of the mass accretion at this stage (see Fig. 1) – Ṁ is smooth and declines after the stars accumulate 95% of the final mass, whereas in the earlier evolution period, Ṁ can be quite variable and feature strong episodic bursts. The effect of variations in this quantity is discussed later in the text. This time instance is marked for every model in Fig. 2 by the red diamonds. The black dots correspond to the evolutionary phase preceding the zero-point age for each object⁴. We refer to this phase as the protostellar phase of evolution. The green symbols mark the reference ages of 1 Myr, 10 Myr, and 25 Myr for each model. The black solid lines present the isochrones for the same ages, but derived from the non-accreting stellar evolution models of (Yorke & Bodenheimer 2008, hereafter, the non-accreting isochrones).

Fig. 3 HR diagrams for models with an age of 1 Myr (top subplot), 10 Myr (center), and 25 Myr (bottom). The black solid lines indicate non-accreting isochrones for 1 Myr (top subplot), 10 Myr (center), and 25 Myr (bottom). The black dashed lines show the non-accreting isomasses. Red and blue markers indicate accreting and non-accreting models, respectively, at the age of 1 Myr (stars), 10 Myr (circles), 25 Myr (squares), and with a mass (in solar masses) shown next to each isomass they fit. Black arrows show the difference between accreting and non-accreting models with the same mass and age.

Evidently, young objects occupy the upper region of the L_∗–T_eff diagram, shifting toward the bottom left part as they age. Notable excursions of young protostellar objects (black dots) to the upper left corner of the diagram are caused by the strong mass accretion bursts discussed in Sect. 2.3. During these bursts, the stellar luminosity and effective temperature both increase. The contraction of the star after the accretion burst occurs on the thermal relaxation timescale t_therm, which can be defined in our case as (16)where the summation is performed over the burst duration defined as a time interval when α_i ≠ 0 (or Ṁ> 10^-5M_⊙ yr^-1). The expression under the sum gives the total accretion energy absorbed by the star during the burst, and this energy is divided by the mean stellar photospheric luminosity averaged over the duration of the relaxation period. We found that the stars recover approximately the previous equilibrium on the thermal relaxation timescales. Our analysis of t_therm shows that the star can spend from several hundreds to up to ten thousand years (depending on the burst strength) in the peculiar excursion tracks, where it can potentially be confused with more massive stars in quiescence. A more detailed study of this phenomenon is deferred to a follow-up paper.

A visual inspection of Fig. 2 indicates that young objects show a notable deviation from the non-accreting isochrones. This is especially evident for the models with effective temperatures log T_eff ≥ 3.5, which lie notably lower than the 1.0 Myr old isochrone, meaning that these objects appear older on the L_∗–T_eff diagram than they truly are⁵. Our results agree with findings of Hosokawa et al. (2011), which state that the non-accreting isochrones can sometimes overestimate stellar ages for stars with effective temperatures above ≥ 3500 K (log T_eff ≈ 3.54). Our results agree in general with numerical models of Baraffe et al. (2016), but their low-mass objects usually show stronger deviations from the non-accreting isochrones.

The non-accreting isochrones are often used to derive the stellar ages using the measured bolometric luminosities and effective temperatures. Our modeling demonstrates that this practice needs to be taken with care, especially for young stars with effective temperatures ≥ 3200 K, which may look older than they truly are when using the non-accreting isochrones. For older stars, however, we obtained a much better fit with the 10 Myr old and 25 Myr old non-accreting isochrones.

It is important to compare our models with the non-accreting isochrones not only for objects of similar age, but also for objects of similar mass. We illustrate this in Fig. 3, where we plot both the isochrones and isomasses derived from the non-accreting stellar evolution models of Yorke & Bodenheimer (2008). In particular, the solid lines in the top, middle, and bottom panels show the non-accreting isochrones for the 1.0 Myr, 10 Myr, and 25 Myr old stars, respectively. The dashed lines in all panels present the isomasses for several chosen values of the stellar mass. The intersection between a specific isochrone and a specific isomass represents a stellar or substellar object with the corresponding age and mass. These intersections are marked with blue symbols. The red symbols denote objects with the same mass and age, but derived from our numerical models. The black arrows indicate the difference between our accreting models and non-accreting models of Yorke & Bodenheimer for similar mass and age. The resulting relative deviations (in per cent) in total luminosities, effective temperatures, and stellar radii of accreting models from the non-accreting ones are summarized in Table 2 for 10 objects shown in Fig. 3. The deviations are calculated using the following formula: (17)where X_accr and X_non−accr stand for L_∗, T_eff, or R_∗ in accreting and non-accreting models, respectively. For the 1.0 Myr old objects, we also calculated the deviations assuming the stellar zero age at the time instance when protostars accumulate 98% of their final mass (in contrast to the originally adopted value of 95%). This new value implies a disk-to-star mass ratio of 2%, which is more appropriate for somewhat more evolved T Tauri stars.

Evidently, notable differences exist (on the order of tens of per cent) between the accreting and non-accreting 1.0 Myr old models in terms of the total luminosity for models of all masses, including the very-low-mass stars and brown dwarfs. This cannot be easily seen from Fig. 2, which in fact suggests a relatively good fit between the accreting and non-accreting models for objects with log T_eff< 3.5. This is a mere coincidence, however, and the accreting and non-accreting models of the same mass are shifted along (but not away from) the non-accreting isochrone. This example illustrates the importance of using both the stellar ages and masses when making a comparison of accreting and non-accreting models, because the accreting objects may fall on the non-accreting isochrone, but be displaced along the isochrone when compared with non-accreting objects of the same mass. All accreting models except for the least massive ones have lower luminosities than their non-accreting counterparts, as can be seen from the sign of the calculated differences in Table 2.

A disagreement of similar magnitude between accreting and non-accreting models is also found for the stellar radius, while the corresponding effective temperatures are characterized by a much smaller deviation. On the other hand, the 10 Myr old and 25 Myr old accreting and non-accreting models show a much better agreement with each other. We note that the accreting models with the 98% zero-age definition show on average higher deviations than the non-accreting models of the same mass and age.

Table 2 provide the differences between accreting and non-accreting models only for ten selected models. In order to check if the found tendency takes place for all models in our sample, we present in Table 3 the values of the differences averaged (by absolute value) over models with M_{∗ ,fin}> 0.2 M_⊙ and models with M_{∗ ,fin}< 0.2 M_⊙. In agreement with our previous conclusions, the largest deviations are found for the total luminosities and stellar radii and the smallest deviations are found for the effective temperatures. In general, the deviations diminish as the stars age, and the averaged differences for the 10 Myr old and 25 Myr old objects are almost an order of magnitude smaller than those for the 1.0 Myr old objects.

Fig. 4 Similar to Fig. 2, but for the hot accretion scenario.

5. Hot accretion

In this section, we present the results for the hot accretion scenario, in which α is set to 0.1 during the entire evolution period. Figure 4 shows the stellar evolutionary sequences for all 35 models. The notations are the same as in Fig. 2. A visual comparison of Figs. 2 and 4 indicates that the hot accretion models behave similarly to the hybrid accretion models. The 1.0 Myr old objects show a moderate mismatch with the corresponding 1.0 Myr old non-accreting isochrone. On the other hand, the 10 Myr old and 25 Myr old objects show a good agreement with the corresponding non-accreting isochrones.

As with the hybrid accretion scenario, we quantify the disagreement between the accreting and non-accreting models and show in Table 4 the relative differences in the total luminosity (ΔL_∗), effective temperature (ΔT_eff), and stellar radius (ΔR_∗) for objects with the same mass and age. For the sake of comparison, we chose the same ten objects as in Table 2. Again, a moderate deviation exists between the accreting and non-accreting 1.0 Myr old models in terms of the total luminosity and stellar radius, but this disagreement diminishes for 10 Myr old and 25 Myr old models. At the same time, the deviation in the effective temperature is rather small for all models, regardless of their mass and ages. All in all, the behavior of the hybrid and hot accretion models is similar, in agreement with the previous work of Baraffe et al. (2012). However, one difference between our work and that of Baraffe et al. is worth emphasizing: for the lowest stellar masses, our accreting models at 1 Myr predict stellar luminosities and radii that are higher than those of non-accreting isochrones (see Tables 2 and 4), while in Baraffe et al. the accreting models yield lower L_∗ and R_∗. We note that for stellar masses ≥ 0.15−0.2 M_⊙, both works agree and predict lower L_∗ and R_∗ for accreting models than for their non-accreting counterparts. We note that the STELLAR non-accreting isochrone shows a sudden turn toward lower luminosities at around T_eff = 3500 K (corresponding to ≈ 0.15−0.2 M_⊙), and this feature is not seen in the Lyon non-accreting isochrones, which can partly explain the mismatch we found. However, the exact reasons require further investigation and comparison between the two stellar evolution codes (STELLAR and Lyon), which is beyond the scope of this work.

Finally, in Table 5 we show the averaged differences in L_∗, T_eff, and R_∗ calculated separately for all models with M_{∗ ,fin}> 0.2 M_⊙ and all models with M_{∗ ,fin}< 0.2 M_⊙. Here, we also find no systematic and significant differences between the hybrid accretion and hot accretion models.

6. Cold accretion

In this section, we present our results for the cold accretion scenario, in which α is always set to a low value of 10^-3 independent of the actual value of the mass accretion rate, meaning that almost all accretion energy is radiated away and only a tiny fraction is absorbed by the protostar. Figure 5 shows the stellar evolutionary sequences for all 35 models of Table 1. The meaning of the colored symbols and lines is the same as in Fig. 2. The disagreement between our accreting models and the non-accreting isochrones is evident not only for the 1.0 Myr old objects, but also for the 10 Myr old objects. In particular, several 1.0 Myr old accreting objects lie near the 10 Myr old non-accreting isochrone. Notable differences with the non-accreting isochrones are even present for the 25 Myr old objects with log T_eff> 3.55.

The cold accretion objects experiencing accretion bursts show strong surges in the luminosity, moving to the upper part of the HR diagram. The reason is that almost all accretion energy is radiated away and the star does not expand dramatically during the burst, having similar effective temperatures and radii during and between the bursts (see Fig. 7). We note that the hybrid or hot accretion objects experience an increase in both luminosity and effective temperature, moving to the upper left part of the HR diagram. The increase in luminosity during the bursts is notably stronger, however, sometimes exceeding 1000 L_⊙, than with hybrid or hot accretion, in which case the stellar luminosity does not exceed 400 L_⊙. This difference can again be explained by the compact size of the cold accretion objects during the bursts, which results in higher accretion luminosities than those of more bloated hybrid or hot accretion objects. According to the recent review by Audard et al. (2014), the strongest luminosity outbursts are observed in FU Ori (340–500 L_⊙), Z Cma (400–600 L_⊙), and V1057 (250–800 L_⊙). None of these object’s luminosities exceeds 1000 L_⊙. We note that our burst models show a good agreement with observations in terms of the model accretion rates (Vorobyov & Basu 2015). Therefore, the lack of very strong luminosity outbursts in FU Orionis-type objects may be an indirect evidence in favor of the hybrid or hot scenario for the thermal efficiency of accretion.

Fig. 5 Similar to Fig. 2, but for the cold accretion scenario.

Fig. 6 Similar to Fig. 3, but for the cold accretion scenario.

In Fig. 6 we compare the total luminosities and effective temperatures of ten selected accreting models with the non-accreting models of Yorke & Bodenheimer for the same mass and age. The meaning of the lines and symbols is the same as in Fig. 3. A comparison of Figs. 3 and 6 shows that the 1.0 Myr old cold accretion models are characterized by much stronger deviations from the corresponding non-accreting models than was found for the hybrid accretion case. More importantly, the cold accretion models show notable deviations at 10 Myr, and for upper-mass objects, even at 25 Myr. To quantify the disagreement between the accreting and non-accreting models, the relative deviations in the total luminosity, effective temperature, and stellar radius for objects with the same mass and age are calculated using Eq. (17) and are shown in Table 6. In addition, Table 7 summarizes the averaged differences in L_∗, T_eff, and R_∗ calculated separately for all models with M_{∗ ,fin}> 0.2 M_⊙ and all models with M_{∗ ,fin}< 0.2 M_⊙.

On average, the disagreement between the cold accretion models and the corresponding non-accreting models of Yorke & Bodenheimer has increased by about a factor of several as compared to the case of hybrid accretion. More specifically, the deviations in L_∗ for the 1.0 Myr old models have grown by a factor of 2, while the deviations in T_eff and R_∗ have grown by factors of 3–4. The disagreement for older models is even more notable: the deviation in L_∗ for the 10 Myr old objects has increased on average by a factor of 4 and may reach as much as a factor of 10 for some upper mass models. The same tendency is found for the stellar radius.

6.1. Comparison of models with different thermal efficiencies of accretion

To better understand the differences in the evolutionary tracks of models with cold, hybrid, and hot accretion, we plot in Fig. 7 the time behavior (from top to bottom) of the total luminosity L_∗, photospheric luminosity L_ph, stellar radius R_∗, and effective temperature T_eff for model 18 without strong accretion bursts (left column) and model 23 with strong accretion bursts (right column). The red, green, and blue lines correspond to cold, hot, and hybrid accretion. The vertical dotted lines mark the zero-point ages when protostars accumulate 95% of their final mass. The stellar radius in the hot or hybrid accretion scenarios (especially in the early evolutionary stages) tends to be significantly greater than that in the cold accretion scenario, in agreement with previous studies (e.g., Baraffe et al. 2012). This trend can be explained by stellar bloating caused by a fraction of accretion energy absorbed by the star in the hybrid and hot accretion scenarios. A similar behavior is evident for the photospheric luminosity, which is higher in the hybrid or hot accretion cases thanks to a larger stellar radius. We note that the mass accretion onto the star changes the properties of the star even in the cold accretion case. The star becomes more compact than in the non-accretion case. Despite the accretion luminosity, the star may still be significantly below the non-accreting isochrone.

Fig. 7 Time behavior (from top to bottom) of the total luminosity, photospheric luminosity, stellar radius, and effective temperature in model 18 without accretion bursts (left column) and model 23 with strong accretion bursts (right column). The red, green, and blue lines correspond to the cold, hot, and hybrid accretion scenarios. The vertical dotted lines define the zero-point ages when protostars accumulate 95% of their final mass.

The behavior of the total luminosity (the sum of the photospheric and accretion luminosities) is somewhat more complex. In the non-burst case, the cold accretion model is initially characterized by a higher total luminosity than their hot and hybrid accretion counterparts, mainly because of a smaller stellar radius, and as a consequence, a higher accretion luminosity in the cold accretion case⁶. In the subsequent evolution, however, the situation changes and the total luminosity in the hot and hybrid accretion models becomes dominant over that in the cold accretion model. The reason is that the accretion luminosity diminishes with time and the photospheric luminosity (which is greater in the hot and hybrid accretion models) starts to provide a major input to the total luminosity in the T Tauri phase (Elbakyan et al. 2016). In the burst case, the total luminosity is also on general is smaller in the cold accretion model than in the hot and hybrid models with one exception – L_∗ is greater in the cold accretion model during the bursts thanks to a smaller R_∗, and as a consequence, higher accretion luminosity. These trends in L_∗ explain why the cold accretion models in Fig. 5 demonstrate greater deviations from the non-accreting isochrones than their hot and hybrid accretion counterparts. The systematically lower total luminosities in models with cold accretion in the late evolution phases result in higher ΔL_∗, as is evident in Tables 6 and 7.

The behavior of the effective temperature is also distinct in the cold accretion and hybrid or hot accretion models. The latter demonstrate sharp increases in T_eff during the bursts, while the former are characterized by a rather smooth evolution of T_eff. As a result, the hybrid or hot accretion objects move to the upper left part of the L_∗–T_eff diagram during the burst, while the cold accretion objects lie in the upper part. This difference in the behavior of the cold vs. hybrid or hot accretion models can in principle be used to constrain the accretion scenario in the protostellar phase of evolution, provided that accretion bursts are sufficiently frequent. We note, however, that the hybrid or hot accretion objects can be confused with more massive stars in quiescence, so that care should be taken to distiuguish the low-mass stars in the burst phase from the higher-mass stars in quiescence. Finally, we note that hot accretion with α ≥ 0.2 fails to explain very low-luminosity objects or VeLLOs (Vorobyov et al. 2017), so that a certain care should be taken when considering hot accretion as a viable scenario.

Fig. 8 Similar to Fig. 2, but for the zero stellar age defined as the time instance when the protostar accumulates 90% of their final mass.

7. Discussion and model caveats

7.1. Postprocessing using pre-calculated accretion rate histories

In our models, we used the precalculated accretion rate histories derived from numerical hydrodynamics simulations of disk formation and evolution (Vorobyov & Basu 2010, 2015). These numerical models were shown to reproduce the main known accretion properties of young stellar objects, for example, the Ṁ–M_∗ steep dependence (Vorobyov & Basu 2008, 2009a) and accretion burst frequency and amplitudes (Vorobyov & Basu 2015). We do not expect that self-consistent simulations coupling the stellar evolution model with the disk hydrodynamics simulations, similar to what has been done in Baraffe et al. (2016), will significantly affect our conclusions because the irradiation by the central star has been taken into account when computing the accretion rates. Nevertheless, we plan these self-consistent studies for the near future.

7.2. The zero stellar age

The effect of variations in the adopted definition for the zero stellar age has been considered in Sect. 4. As Table 2 demonstrates, assuming the stellar zero age at the time instance when protostars accumulate 98% of their final mass (in contrast to the value of 95% used throughout the paper) results in a stronger disagreement between the accreting models and the non-accreting isochrones. In the opposite case of a smaller fraction of the accreted mass, we may expect a better agreement. Figure 8 shows the model tracks in the hybrid accretion model, but for the zero stellar age (the red diamonds) defined as the time instance when the protostar accumulate 90% of their final mass. This definition of the zero age corresponds to the evolutionary phase when 10% of the final stellar mass is still confined in the disk and residual envelope, which is more appropriate to the late Class I phase of protostellar evolution than to the T Tauri phase. As Fig. 1 shows, accretion rates at these times may still show some substantial variability. A visual comparison of Figs. 2 and 8 demonstrates that the agreement of the 1.0 Myr old accreting objects (the green stars) with the corresponding non-accreting isochrone is now better than in the case 95% case. However, the zero-age (or the birth) locations are now widely scattered, because the protostars continue accreting at this stage and may experience significant excursions in the value of L_∗ because of the time-varying accretion rates in the late protostellar phase. Time variations in the birth locations of stars present no problem and may be a real phenomenon, as was noted in Baraffe et al. (2012). The classic smooth birthline of stars introduced in Stahler (1983) was calculated based on spherically symmetric collapse simulations adopting a constant protostellar accretion rate and neglecting the possible significant effect of accretion variability introduced by the circumstellar disk⁷. However, the 90% definition may still technically correspond to the late protostellar phase (rather than T Tauri phase), and taking the 90% case as the zero-age point is therefore less well physically motivated than the other two cases.

7.3. The thermal efficiency of accretion

In this work, we have adopted α = 0.1 for the maximum fraction of absorbed energy by the central object and Ṁ_cr = 10^-5M_⊙ yr^-1 for the transition from cold to hot accretion (see Sect. 3.2). The possible variations in these quantities may somewhat affect the stellar tracks on the L_∗–T_eff diagram. For instance, by increasing the value of Ṁ_cr, we simultaneously decrease the time during which the star may absorb part of the accretion energy. In the limit of Ṁ_cr ≫ 10^-4M_⊙ yr^-1, the hybrid accretion tracks will behave similar to those of the cold accretion scenario because accretion bursts of such high magnitude are extremely rare. In the opposite case of Ṁ_cr ≪ 10^-6M_⊙ yr^-1, the hybrid accretion track will resemble those of the hot accretion scenario because the star will spend most of its time absorbing part of the accretion energy. In addition, varying the value of α may produce a spread in the L_∗–T_eff diagram, as was shown in Baraffe et al. (2012). In the hot accretion scenario, however, values of α greater than 0.2 are unlikely because the resulting photospheric luminosities are always greater than 0.1 L_⊙. This contradicts the apparent existence of very low-luminosity objects with the internal luminosity < 0.1 L_⊙ in the protostellar phase (see Vorobyov et al. 2017).

7.4. Initial stellar radius and mass

To check the effect of the initial conditions imposed on the protostellar seed, we recalculated some typical models in the hot and cold accretion scenarios by varying the initial seed radius from 2.7 to 4.0 Jovian radii and reducing the initial seed mass to 4.0 Jovian masses. The resulting pre-main-sequence evolution showed no significant deviation from the evolution of the original models. Unfortunately, we were not able to vary the initial seed parameters over a wider range of radii and masses because of numerical divergence of the stellar evolution code. We plan to perform a more rigorous study of the effects of initial conditions in future self-consistent simulations.

8. Conclusions

We have considered the pre-main-sequence evolution of low-mass stars and brown dwarfs starting from the formation of a protostellar seed and taking into account the mass accretion during the initial several Myr of evolution. The stellar evolution was computed using the STELLAR evolution code originally developed by Yorke & Bodenheimer (2008) and further modified by Hosokawa et al. (2013).

The mass accretion rates were taken from numerical hydrodynamics simulations of disk evolution starting from the gravitational collapse of prestellar cloud cores of various mass and angular momentum. The resulting accretion rates exhibit various patterns of time behavior: from smoothly declining rates to strongly time-varying rates featuring episodic bursts similar in magnitude to those of the FU Orionis-type eruptions (Vorobyov & Basu 2010, 2015). Three scenarios for the thermal efficiency of accretion were considered: hot accretion with a constant fraction of accretion energy absorbed by the central object (star or brown dwarf), cold accretion in which essentially all accreted energy is radiated away, and hybrid accretion in which the fraction of absorbed energy depends on the accretion rate. We compared the resulting stellar evolution tracks in the total luminosity (L_∗) – effective temperature (T_eff) diagram with the isochrones and isomasses derived using the non-accreting models of Yorke & Bodenheimer. Our key findings can be summarized as follows:

• In the hybrid accretion case, young 1.0 Myr old ob-jects show notable deviations from the non-accreting isochronesand isomasses for both low-mass stars and brown dwarfs. Thelargest deviations (relative to the non-accreting models) arefound for the total stellar luminosity and stellar radius (19–51%and 8–32%, respectively), while for the effective temperature thedeviations are relatively mild (0.15–2.1%). The disagreementbetween the accreting and non-accreting models diminisheswith stellar age, remaining within several percent for 10 Myr old and 25 Myrold objects. The calculated deviations depend somewhat on ourdefinition of the stellar zero age.

• The hot accretion case is qualitatively similar to hybrid accretion, but showing somewhat smaller deviations for L_∗ and R_∗.

• The cold accretion case features the largest deviations from the non-accreting models of Yorke & Bodenheimer. For the 1.0 Myr old objects, the deviations in L_∗ in R_∗ are on average factors of 2–3 greater than in the cases of hybrid or hot accretion. A disagreement is also found for the 10 Myr old objects and even for the 25 Myr old objects, especially for the upper mass objects. The greater disagreement as compared to the hybrid or hot accretion cases can be explained by the systematically smaller stellar radii, and as consequence, photospheric luminosities (which dominate the accretion luminosity in the T Tauri phase, Elbakyan et al. 2016) in models with cold accretion.

• As a result of this mismatch, the use of the L_∗–T_eff diagram may lead to the false age estimate for objects with T_eff> 3500 K, as was also previously noted in Baraffe et al. (2009) and Hosokawa et al. (2011). For instance, 1.0 Myr old objects with cold accretion can be falsely identified as 4.5 Myr old ones. For the case of hybrid and hot accretion, the error remains within a factor of 2 for object of 1.0 Myr age, but diminishes for older objects.

• Hybrid and hot accretion models show sharp increases in both L_∗ and T_eff during accretion bursts. As a result, these models show notable excursions to the upper left region of the L_∗–T_eff diagram. In contrast, T_eff in the cold accretion scenario is weakly affected by the bursts, and the cold accretion objects show strong surges only in L_∗. These differences between the stellar evolution tracks of the hybrid or hot accretion and cold accretion objects experiencing accretion bursts can potentially be used to constrain the thermal efficiency of accretion, but care should be taken to avoid confusing the low-mass stars in the burst phase with the higher-mass stars in quiescence.

• The increase in luminosity of the cold accretion objects experiencing accretion bursts is notably stronger than that of the hybrid or hot accretion objects. In particular, the cold accretion objects sometimes exceed 1000 L_⊙ in luminosity, whereas the hybrid or hot accretion objects never exceed 400 L_⊙ in luminosity. According to the recent review by Audard et al. (2014), the strongest luminosity outbursts are observed in FU Ori (340–500 L_⊙), Z CMa (400–600 L_⊙), and V1057 (250–800 L_⊙), all of which are below 1000 L_⊙. Therefore, the lack of very strong luminosity outbursts in the FU Orionis-type objects may be an indirect evidence in favor of the hybrid or hot scenarios for the thermal efficiency of accretion. The hot scenario with α ≥ 0.2, however, fails to reproduce the very low-luminosity objects (VeLLOs) because the resulting photospheric luminosity (which sets the floor for the total luminosity in the early stages of evolution) is always greater than the VeLLO limit of 0.1 L_⊙ (Vorobyov et al. 2017).

In the appendix, we provide the updated isochrones based on the stellar evolution models taking mass accretion into account. The usefulness of these isochrones hinges on the adopted scenario for the thermal efficiency of accretion, and efforts now should be made to find the means of distinguishing between the hybrid or hot and cold accretion cases, similar to what has recently been reported in Geroux et al. (2016).

Acknowledgments

We are thankful to the anonymous referee for helpful suggestions that allowed us to improve the manuscript. E.I.V. and V.G.E. are thankful to Isabelle Baraffe and Gilles Chabrier for stimulating discussions that inspired this work. E. I. Vorobyov and V. G. Elbakyan acknowledge support by the Russian Ministry of Education and Science Grant 3.5602.2017. V. G. Elbakyan acknowledges the Southern Federal University for financial support with an international travel grant. The simulations were performed on the Vienna Scientific Cluster (VSC-2 and VSC-3) and on the Shared Hierarchical Academic Research Computing Network (SHARCNET). We also appreciate the financial support by the Grants-in-Aid for Basic Research by the Ministry of Education, Science and Culture of Japan (16H05996: TH) and by Grant-in-Aid for JSPS Fellows (SH). Portions of this work were conducted at the Jet Propulsion Laboratory, California Institute of Technology, operating under a contract with the National Aeronautics and Space Administration (NASA).

References

Appendix A: Updated isochrones using accreting models

Fig. A.1 Total luminosity–effective temperature diagram at 1.0 Myr and 10 Myr in the hybrid, hot, and cold accretion scenarios. The colored symbols represent our accreting models and the color of the symbols varies according to the stellar mass shown in the vertical bar (in M⊙). The corresponding ages of the model data and scenarios are indicated in each panel. The black dashed lines are the isochrones (of the corresponding age) derived from the non-accreting models of Yorke & Bodenheimer (2008), while the black solid lines are the best-fit curves to our model data (accreting isochrones, see the text). The red solid lines represent the non-accreting isochrones of Yorke & Bodenheimer that fit our model data best.

In this section, we calculate the isochrones taking accretion into account and compare them with the isochrones derived from the non-accreting models of Yorke & Bodenheimer (2008). The colored symbols in Fig. A.1 represent our accreting models in the L_∗–T_eff diagram at 1.0 Myr and 10 Myr in the hybrid, hot, and cold accretion scenarios. The corresponding ages and scenarios are indicated in the panels. The color of the symbols varies according to the stellar mass shown in the vertical bar (in M_⊙). The black dashed lines are the corresponding isochrones derived from the non-accreting models of Yorke & Bodenheimer (2008), while the black solid lines are the best-fit curves to our model data (hereafter, accreting isochrones). The interpolation is made using the Vandermonde matrix and a forth-order polynomial of the form (A.1)and the coefficients of the polynomical are listed in Table A.1.

Evidently, a notable mismatch exists between the non-accreting (black dashed lines) and accreting (black solid lines) isochrones, especially for the cold accretion objects of 1 Myr age. To quantify this disagreement, we plot in Fig. A.1 with the red solid lines the non-accreting isochrones that fit our model data best. The resulting ages of these isochrones are 1.5 Myr (top left panel), 1.5 Myr (top right panel), 4.5 Myr (bottom left panel), and 16 Myr (bottom right panel). As was discussed earlier, the largest error in the age determination when using the non-accreting isochrones can occur for young, intermediate, and higher-mass objects with cold accretion – the 1.0 Myr old objects can be falsely identified as 4.5 Myr old objects. The error still remains significant at 10 Myr – these objects can be falsely interpreted as 16 Myr old objects. For the case of hybrid and hot accretion, the error remains within a factor of 2 for object of 1.0 Myr age, but diminishes for older objects.

In order for the accreting isochrones to be of practical use, efforts should now be focused on determining the most plausible scenario for the thermal efficiency of accretion. Multidimensional numerical hydrodynamics simulations of accreting stars, similar to what has recently been reported by Geroux et al. (2016), can give us a clue about the thermal efficiency of accretion, but these simulations are usually too computationally expensive to be used to simulate the long-term accretion history of (sub-)solar mass stars. Nevertheless, they can be useful in calibrating the less computationally expensive one-dimensional models of stellar evolution.

All Tables

All Figures

Fig. 1 Mass accretion rates vs. time for 35 models. The horizontal dashed lines mark the critical value for the transition from cold to hot accretion scenario. The models are arranged so that the final stellar mass increases from left to right and from top to bottom. The vertical dashed, dotted, and dash-dotted lines mark time instances when 90%, 95%, and 98% of the total stellar mass has been accumulated.

In the text

Fig. 2 Stellar evolution sequences on the total luminosity L∗ – effective temperature Teff diagram for the hybrid accretion models. The dots present the model tracks and the color of the dots varies according to the stellar age shown in the vertical bar (in log yr). The zero-point age for each model is marked by the red diamonds. The green symbols mark the reference ages, as indicated in the bottom left corner, that has elapsed since the zero-point age of each object. The black solid lines provide the isochrones for stellar ages of 1 Myr, 10 Myr, and 25 Myr (from top to bottom) derived from the non-accreting stellar evolution models of Yorke & Bodenheimer (2008).

In the text

Fig. 3 HR diagrams for models with an age of 1 Myr (top subplot), 10 Myr (center), and 25 Myr (bottom). The black solid lines indicate non-accreting isochrones for 1 Myr (top subplot), 10 Myr (center), and 25 Myr (bottom). The black dashed lines show the non-accreting isomasses. Red and blue markers indicate accreting and non-accreting models, respectively, at the age of 1 Myr (stars), 10 Myr (circles), 25 Myr (squares), and with a mass (in solar masses) shown next to each isomass they fit. Black arrows show the difference between accreting and non-accreting models with the same mass and age.

In the text

	Fig. 4 Similar to Fig. 2, but for the hot accretion scenario.
In the text

	Fig. 5 Similar to Fig. 2, but for the cold accretion scenario.
In the text

	Fig. 6 Similar to Fig. 3, but for the cold accretion scenario.
In the text

Fig. 7 Time behavior (from top to bottom) of the total luminosity, photospheric luminosity, stellar radius, and effective temperature in model 18 without accretion bursts (left column) and model 23 with strong accretion bursts (right column). The red, green, and blue lines correspond to the cold, hot, and hybrid accretion scenarios. The vertical dotted lines define the zero-point ages when protostars accumulate 95% of their final mass.

In the text

	Fig. 8 Similar to Fig. 2, but for the zero stellar age defined as the time instance when the protostar accumulates 90% of their final mass.
In the text

Fig. A.1 Total luminosity–effective temperature diagram at 1.0 Myr and 10 Myr in the hybrid, hot, and cold accretion scenarios. The colored symbols represent our accreting models and the color of the symbols varies according to the stellar mass shown in the vertical bar (in M⊙). The corresponding ages of the model data and scenarios are indicated in each panel. The black dashed lines are the isochrones (of the corresponding age) derived from the non-accreting models of Yorke & Bodenheimer (2008), while the black solid lines are the best-fit curves to our model data (accreting isochrones, see the text). The red solid lines represent the non-accreting isochrones of Yorke & Bodenheimer that fit our model data best.

In the text