Explaining the luminosity spread in young clusters: proto and pre-main sequence stellar evolution in a molecular cloud environment
MMNRAS , 1–20 (2017) Preprint 1 February 2018 Compiled using MNRAS L A TEX style file v3.0
Explaining the luminosity spread in young clusters:proto and pre-main sequence stellar evolution in amolecular cloud environment
Sigurd S. Jensen (cid:63) & Troels Haugbølle † Centre for Star and Planet Formation, Niels Bohr Institute and Natural History Museum of Denmark, University of Copenhagen,Øster Voldgade 5-7, DK-1350 Copenhagen K, Denmark
Accepted 2017 October 31. Received October 30; in original form 2017 June 27
ABSTRACT
Hertzsprung-Russell diagrams of star forming regions show a large luminosity spread.This is incompatible with well-defined isochrones based on classic non-accreting pro-tostellar evolution models. Protostars do not evolve in isolation of their environment,but grow through accretion of gas. In addition, while an age can be defined for a starforming region, the ages of individual stars in the region will vary. We show how thecombined effect of a protostellar age spread, a consequence of sustained star formationin the molecular cloud, and time-varying protostellar accretion for individual proto-stars can explain the observed luminosity spread. We use a global MHD simulationincluding a sub-scale sink particle model of a star forming region to follow the accre-tion process of each star. The accretion profiles are used to compute stellar evolutionmodels for each star, incorporating a model of how the accretion energy is distributedto the disk, radiated away at the accretion shock, or incorporated into the outer lay-ers of the protostar. Using a modelled cluster age of 5 Myr we naturally reproducethe luminosity spread and find good agreement with observations of the Collinder 69cluster, and the Orion Nebular Cluster. It is shown how stars in binary and multiplesystems can be externally forced creating recurrent episodic accretion events. We findthat in a realistic global molecular cloud model massive stars build up mass over rel-atively long time-scales. This leads to an important conceptual change compared tothe classic picture of non-accreting stellar evolution segmented in to low-mass Hayashitracks and high-mass Henyey tracks.
Key words: stars: formation, stars: luminosity function, stars: pre-main-sequence,open clusters and associations: individual: Collinder 69, Orion Nebular Cluster
Traditionally, the stellar structure of new-born stars havebeen calculated as a closed system, assuming initially a lowsurface temperature, large radii, and a given mass. Placedon the Hayashi track the protostar first contracts at almostconstant effective temperature, and then moves towards themain sequence. In reality, stars are born from collapsingcores in molecular clouds and grow through accretion of gason to a circumstellar accretion disk, where matter is sub-sequently transported to the star. The accretion rate startsout at a relatively high pace and slows down while mat-ter in first the envelope, during the protostellar phase, andthen in the circumstellar disk, during the pre-main sequence (cid:63)
Contact e-mail: [email protected] † Contact e-mail: [email protected] phase (PMS), is exhausted. The accretion is not smoothlydecreasing but rather intermittent and unsteady in nature.Observations of embedded protostars as well as FU Orionisand Ex-Lupin type PMS stars show convincing evidence ofdynamical accretion both during the protostellar and PMSphase (Herbig 1966, 1977, 1989; Kenyon et al. 1990; K´osp´alet al. 2007; Evans et al. 2009; Scholz et al. 2013; Jørgensenet al. 2015; Safron et al. 2015; Hunter et al. 2017). Further-more, young clusters are observed to have a large spread inthe observed luminosity in the Hertzsprung-Russell (HR) di-agram. Classical non-accreting formation models lead to welldefined isochrones, which have been used for decades for ageestimation of star forming regions. However, the large lumi-nosity spread observed in such populations has been a causefor concern. In the classical paradigm such large spreadsin luminosity imply active star formation over tens of mil-lions of years and thus favour a slow star formation process, © a r X i v : . [ a s t r o - ph . S R ] J a n Jensen & Haugbølle which is in conflict with observed and modelled lifetimes ofstar formation regions (Padoan et al. 2016). Baraffe et al.(2009, 2012) suggested that the apparent age spread in starforming regions could be explained by accreting star formingmodels, specifically models featuring accretion bursts.A central feature of accreting models is the energy re-leased at the accretion shock in the boundary layers of theprotostar. Potential energy is heating the disk through vi-cious processes and becomes kinetic energy when matter istransported towards the star. At the accretion shock the in-fall is stalled and kinetic energy turns into thermal energy,which is either radiated away or absorbed by the protostar.The first scenario is termed cold accretion while the sec-ond scenario where a fraction of the accretion luminosityheats the protostar is known as warm or hot accretion. Whatfraction of the accretion luminosity is absorbed by the starremains an open question and probably depends on the con-ditions of the shock front and the geometry of the accretionflow. Hot and cold accretion models behave very differentlyand robust models must feature a reasonable approximationfor the energy dynamics at the shock front.We present a new framework to understand the lumi-nosity spread in young stellar clusters. As input for the stel-lar evolution models we use a realistic large-scale simulationof a molecular cloud fragment (4 pc) in size that contains alarge sample of 413 protostars, with precisely recorded massaccretion rates. We evolve each protostar from the birth as asecond Larson core to the main sequence. To investigate theimpact of the detailed accretion process we run a number ofdifferent stellar structure evolution models for each proto-star, which include different degrees of cold and hot accre-tion, including an accretion rate dependent model. The HRdiagrams for the synthetic clusters are compared to obser-vations of the Collinder 69 and Orion Nebular Cluster starforming regions, and it is shown that a combination of twoeffects serve to explain the observed luminosity spread. Boththe non-steady accretion rate of individual protostars, andthe age spread between cluster members serves to broadenthe HR diagram, compared to classical models, and bring itin agreement with observations. The simulation is carried out with the ramses code (Teyssier2002) modified to include random turbulence driving, sinkparticles as a sub-grid model for protostars, technical im-provements allowing for efficient scaling to several thousandcores, and an improved HLLD solver that is stable in su-personic flows with high Mach numbers. We refer the readerto Haugbølle et al. (2017) for an in-depth description of thesimulation, and the resulting molecular cloud structure. Inthat article, there is also a full account of the sink particlemodel, implemented physics, and code modifications com-pared to the public version of ramses .Our model setup is a periodic box of size (4 pc) witha total mass of 3,000 M (cid:12) and an isothermal temperatureof 10 K. Assuming a molecular weight of 2.37 the averagenumber density in the box is 795 cm − . The initial magneticfield strength is 7.2 µ G. In our model we solve the isother-mal MHD equations. Starting from a uniform gas density and magnetic field, using a solenoidal driving on the largestscales, the RMS velocities are amplified to a sonic Machnumber M = σ v,3D / c s of approximately 10, where σ v, 3D isthe volume averaged 3D RMS velocity and c s is the isother-mal sound speed. To let the turbulent flow develop and reacha statistically steady state we first evolve the box for 20 dy-namical or turn-over times, where t dyn = L box /( σ v,3D ) , be-fore self-gravity and the sink particle recipe is turned on. Theturbulent driving mimics the cascading external forcing fromlarger scales not covered in the model. The average density,RMS velocity, and magnetic field strength in the model aretypical for parsec sized molecular clouds (Heyer et al. 2009;Crutcher 2012). The balance between kinetic and gravita-tional energy, the virial parameter, determines the rate atwhich stars are formed (Padoan et al. 2012, 2017). A typicaldefinition is to use the equivalent virial number for a homo-geneous sphere, α vir = σ v,3D R /( GM ) . For our model thisyields α vir = . , where we have used R = L box / . This isvery similar to observed values for molecular clouds that arein the range 0.5 to 5 (Kauffmann et al. 2013). The relativelylow value of α vir is probably the cause for the late-time for-mation of a gravitationally bound structure in the box.The only difference between the simulation used inthis paper and in previous articles (Padoan et al. 2014b;Frimann et al. 2016b) is the imposition of a small maximumtimestep at the highest level of refinement of 80 days. Atevery timestep we write out the sink particle properties, re-sulting in 150 GB data for 413 stars and a high-resolutionaccount of the accretion rate for each sink particle, with upto 1 million records for the oldest sink particle in the model.The small timestep also facilitates an accurate longtime in-tegration of sink particle orbits in the simulation. Finally,compared to earlier simulations, we have corrected a bug inthe innermost part of the smoothed gravitational potentialused for each sink particle that could in some circumstanceslead to the disruption of systems with bound sinks, whentwo sinks are orbiting inside a single cell.The simulation uses an adaptive mesh with a root gridof 256 and 6 levels of refinement reaching a maximum res-olution of ∆ x min = pc / = AU. This is enough tocapture the formation process of individual stars, and thedynamical creation of binary and multiple system throughgravitational interaction. It is insufficient to resolve circum-stellar disks, and this is important to keep in mind wheninterpreting accretion rates. Accretion rates as reported andused in the paper are in reality infall rates from the sur-rounding gas reservoir to the circumstellar disk. The realintermittency in the accretion rates could either be higher,if disk mediated instabilities transport angular momentumin an uneven fashion from 50 AU down to the star, or lower,if the disk acts as a buffer that smooths out sudden bursts ofinfall as can happen (periodically) in e.g. binary systems. Inany case, infall rates on larger scales supplies the mass thatis accreted on smaller scales, and the two must necessarilyagree in a time-averaged sense. To account for the effects ofthe unresolved outflows we only accrete 50% of the mass andmomentum to the sink particle (see Haugbølle et al. 2017,for an in-depth description of the sink accretion procedure).Refinement is based solely on density, because our maingoal is to resolve the gravitational collapse. We refine anycell at the root grid level that has a density 10 times abovethe average density. Then we add new levels every time the
MNRAS000
MNRAS000 , 1–20 (2017) he luminosity spread in young clusters density increases by a factor of four, in order to keep theminimum number of cells per Jeans length L J = (cid:115) c s π G ρ ∆ x (1)constant. In our case we have L J > . everywhere exceptat the highest level of refinement. A sink particle is createdwhen the gas reaches a density, at the highest level of refine-ment, such that L J < corresponding to a number densityof n s > . cm − . For details of the sink particle recipesee Padoan et al. (2014b).We run the model for 2.55 Myr, from the time whenthe first star is formed. The model has a global free-falltime of t ff = . Myr and a dynamical time t dyn = . .The model has therefore been evolved for a bit more thantwo free-fall or dynamical times with star formation turnedon. This allows for a well evolved stellar population withmany T-Tauri stars, and makes the model comparable inage to some of the youngest nearby star forming regions. Atthe end of the run 214 stars are more than a million yearsold, while 328 stars have an accretion rate that is less than − M (cid:12) yr − , indicating that they have finished their mainaccretion phase. Below we make a more detailed selectioncriteria and find that 182 of the stars are deeply embedded,and therefore in the Class 0 phase. The embedded stars aremostly young and small. 40 have a mass below the browndwarf limit of 0.076 M (cid:12) , and 127 have a mass below 0.3 M (cid:12) . Of the embedded stars, only 16 are more massive than1 M (cid:12) . The median age of the embedded stars is 480 kyr.The star formation efficiency (SFE), star formation rateper free-fall time (SFR ff ) and initial mass function (IMF) areimportant indicators of the numerical fidelity and realism ofthe model. The star formation efficiency is the fraction ofmass turned into stars: SFE = M sink / M box , while the starformation rate per free-fall time is the SFE per free-fall time:SFR ff = dSFE / d ( t / t ff ) . In Fig. 3 is shown the star formationefficiency and star formation rate per free-fall time. Theyare well fitted by a power-law with SFR ff ≈ . ( t / t ff ) . slowly increasing from 0.04 at 0.5 Myr (0,4 t ff ) to 0.08 at 2.4Myr (2 t ff ). This increase in the SFR ff is probably relatedto the formation of a single dense cluster at the end of therun (see Fig. 1 for column densities at different evolutionstages). The increasing density will decrease the local virialnumber increasing the star formation rate (Padoan et al.2012, 2017). Furthermore, at the end of the run the centralcluster is similar in size to the box and therefore the drivingscale, making it impossible to disrupt it further. To continuefor a longer time period, keeping the evolution realistic, ex-ternal forcing with scales beyond the 4 pc box-size wouldhave been needed to interact properly with the cluster, suchas in Padoan et al. (2017). A SFR ff of 4% to 8% is realisticfor the star forming regions considered below (Padoan et al.2014a). In Fig. 2 is shown the IMF sampled at 1.6, 1.9 and2.2 Myr, and an additional IMF extracted from an identicalrun, where the only difference is an 8 times lower overall res-olution. This additional run was used to validate the extentof the numerical convergence of the model. The IMF is ingood agreement with the “Chabrier IMF” (Chabrier 2003),and the low- and intermediate-mass stellar populations areessentially time-independent at late times, with only a verti-cal shift due to the growing total stellar mass and the higherend of the Salpeter slope still being developed. Comparing the IMF at the two different resolutions we estimate the stel-lar population to be complete above (cid:38) . (cid:12) . To resolve thebrown dwarf population would require higher resolution tobetter sample the turbulence at all scales and recover therare very high density peaks, required for the formation ofbrown dwarf mass stars, but for the present purpose thecost of such a run – ∼
10 million cpu hours – is prohibitive.Furthermore, in this paper we are mainly comparing to ob-servations of stars with masses above (cid:38) . (cid:12) , and resolvingthe brown dwarfs population while interesting in itself wouldnot add additional constraints.In this work it is important to distinguish between thecluster age, which is defined as the time since the formationof the first star in the molecular cloud simulation, and thestellar age, which is defined from the creation of the sinkparticle.Below we compare the sink particle population withclusters of an estimated age of ∼ ff is not constantthroughout the simulation, but the time variation in theSFR ff is acceptable. The current run has already reached anSFE of 13%. Stitching two copies of the stellar evolutionaryprofiles together, with the the older profile only containingnon-accreting stars, carries the implicit assumption that starformation continues for up to 5 Myr. A larger box-size wouldbe required to properly explore this, but continuing runningthe model to 5 Myr is, apart from costly – it would cost ∼ ∼ − t dyn = − Myr, and external feed-back and infall from larger scales would play an importantrole in the further evolution of the cloud (Dobbs et al. 2014;Padoan et al. 2016).In the rest of the paper we will refer to sink particles as“stars”. mesa
The evolution of individual protostars are modelled in theone-dimensional stellar evolution code mesa version 8845.We refer to the instrument papers of mesa (Paxton et al.2011, 2013, 2015) for a complete description of the stellarevolution code. The code has been augmented with modulesfor accretion of mass with a variable entropy and diagnosticsto calculate the accretion luminosity. The implementation of
MNRAS , 1–20 (2017)
Jensen & Haugbølle
Figure 1.
Column density from left to right and top to bottom at t = 1.6, 1.9, 2.2 and 2.5 Myr after the formation of the first star.The red circles mark the positions of stars with masses M < . (cid:12) , brown triangles are stars with . (cid:12) < M < . (cid:12) , while orangesquares indicate stars with M > . (cid:12) . Numbers refer to the stars shown in Fig. 4. the accretion luminosity is described in detail below. We donot include rotation in the models.Each protostar is evolved from the same initial setupand with the same physics implementation. The physicalparameters were chosen to be suitable for stars of low- andintermediate masses.Convection is implemented using the mixing length the-ory adopted from Cox & Giuli (1968) with α M LT = . . Weadopted overshooting parameters for low- and intermediate-mass stars from Choi et al. (2016) above and below convec-tive zones. Convective stability is determined based on theSchwarzschild criterion. The Ledoux criterion for convection was tested and no notable difference between the two criteriawere found for the protostellar models. Similarly a numberof alternative implementations of the mixing length theoryincluded in mesa have been tested and found to producesimilar results.The reaction networks for H , Li , Be and B wereenabled in addition to the standard nuclear reaction networkof mesa . mesa allows for different opacity tables in different tem-perature regimes. In the low temperature regime ( T < , K) we use opacity tables from Freedman et al. (2008) and forhigher temperatures the tables of Grevesse & Sauval (1998)
MNRAS000
MNRAS000 , 1–20 (2017) he luminosity spread in young clusters Figure 2.
Initial mass function sampled at t = 1.6, 1.9, and 2.2Myr, corresponding to an SFE of 6%, 8%, and 10%, and the IMFat an SFE of 10% for a simulation with 8 times lower resolu-tion. The histogram bins have been slightly offset for clarity. Thedashed line is a Chabrier IMF normalised to the same total stellarmass as at 2.2 Myr. The vertical shifts in the histograms is due tothe increased stellar mass as a function of time. The only signif-icant difference between the different IMFs is the growth in theSalpeter range with time, which happens because the formationtime of most massive stars is comparable or larger than the ageof the cluster. ff ]0.000.020.040.060.080.100.12 S F E = M s i n k / M bo x ff ]0.000.020.040.060.080.100.12 S F R ff Figure 3.
The evolution of the star formation efficiency and starformation rate per free-fall time (inset) as a function of time, mea-sured in free-fall times, in our model. The dashed line is a power-law fit SFE = . ( t / t ff ) . corresponding to SFR ff = . ( t / t ff ) . . is used. We use the standard mesa equation of state (EOS),which is based on the OPAL EOS tables of Rogers & Nay-fonov (2002) and extended to lower temperatures with thetables of Saumon et al. (1995). The initial setup for the mesa models is a relaxed sec-ond Larson core of uniform composition. Recent studies ofthe collapse of prestellar cores with a broad spectrum of different initial conditions found second Larson cores with M ∼ . (cid:12) and R ∼ . (cid:12) to be universal for low- andintermediate-mass stars (Vaytet & Haugbølle 2017) and wehave chosen to use cores with this mass for our simulations.We have not fixed the initial radius and due to the differ-ences between the opacity tables we have utilised in mesa and those of Vaytet & Haugbølle (2017) we get an initialradius of R = . (cid:12) .The composition of the core is uniform and the massfractions for the essential elements of the protostar are basedon the values for the local interstellar medium (Ferri`ere 2001;Tosi 2000; Prodanovi´c et al. 2010): X = . , Y = . , H = × − and He = . × − . For the heavier elementswe use the mass fractions of Grevesse & Sauval (1998). Weaccrete matter with the same composition as the initial core. The accretion shock near the surface of the protostar is acrucial part of accreting star formation models, but our un-derstanding of the energy transfer at the shock remains lim-ited.We have adopted the terminology of Baraffe et al.(2009) for the energy transfer at the accretion shock. Theaccretion luminosity is divided into two terms, one for theenergy radiated away from the protostar and one for theenergy that is injected into the protostar. L outacc = (cid:15) ( − α ) GM (cid:219) MR L inacc = (cid:15)α GM (cid:219) MR (2)Where M and R are the protostellar mass and radius re-spectively, α is the thermal efficiency parameter, G is thegravitational constant and (cid:219) M is the accretion rate. The fac-tor (cid:15) depends on the details of the accretion process, with (cid:15) (cid:54) for gravitationally bound material and (cid:15) (cid:54) . forboundary layer accretion from a thin disk. We set (cid:15) = . aswe assume thin disk accretion as previously done by Baraffeet al. (2009). This requires that half of the accretion lumi-nosity is stored in rotational energy or lost during accretionin the disk. This would be the case if e.g. the distance be-tween the inner gap of the accretion disk and the stellarsurface is one stellar radius and that it coincides with theco-rotation radius (Hartmann et al. 2016; Hartmann 2008).The size of the gap can be measured indirectly by assum-ing co-rotation between the foot points of magnetosphericaccretion channels originating close to the inner edge andthe PMS star, and support such a size range (Carr 2007;Eisner et al. 2010). In the earliest stages of star formationbefore an accretion disk has formed and has become thedominant means of transferring mass from the envelope tothe protostar this approximation is not valid, which causesthe accretion luminosity we calculate in that case to be alower limit. This is also the case if the foot points of the ac-cretion channels or the co-rotation radius are further awaythan two stellar radii.Spherical accretion without rotation was studied exten-sively by Stahler et al. (1980a,b, 1981) and they found thatfor an optically thick accretion shock the maximum thermalefficiency is 75% . Observations of T-Tauri stars have indi- In the case of (cid:15) = . then α = . MNRAS , 1–20 (2017)
Jensen & Haugbølle cated that a large fraction of the gas accreted in this evolu-tionary stage happens through an accretion disk. Hartmannet al. (1997) argue that accretion through a disk must becold or have a very low thermal efficiency. This is due tothe optically thick nature of the protostar compared to thesurrounding regions, into which the accretion luminosity canradiate freely when accretion occurs on a limited region ofthe protostellar surface through e.g. magnetospheric accre-tion channels.Geroux et al. (2016) recently compared two-dimensionalsimulations of accreting young stars with previous one-dimensional models of accreting protostars. They found areasonable agreement between the different models whilethere was a tendency for the one-dimensional models to ex-aggerate the effects of a high thermal efficiency relative tothe two-dimensional models.Below, we present results for a number of constant ther-mal efficiencies as well as an dynamic α -model in which thethermal efficiency is a function of the instantaneous accre-tion rate. The model is motivated by the assumption thatthe accretion flow must cover the majority of the protostellarsurface during high accretion rates leading to a high ther-mal efficiency, while accretion at lower rates is expected tobe more localised occurring for example as magnetosphericaccretion through a magnetised loop in which case materialwill be aggregated to the stellar surface at a low thermalefficiency with the same entropy as surface material.As a result of the work presented by Geroux et al. (2016)we have limited α (cid:48) . . Together with the choice of (cid:15) = . this means that at most 25% of the accretion energy is ab-sorbed by the protostar. This is a conservative upper limitfor the thermal efficiencies and the effect of hot accretioncan be considered moderate in this case. A higher thermalefficiency leads to a larger protostar as the continuous addi-tion of energy alters the thermal structure of the protostar.This results in a longer contraction phase towards the mainsequence once accretion is stopped and delays the onset ofnuclear processes such as deuterium and lithium burning.The dynamic α ( (cid:219) M ) -model is constructed as a step func-tion with a smooth transition. It is given by the followingexpression α ( (cid:219) M ) = α L exp ( (cid:219) M m ∆ ) + α H exp ( (cid:219) M ∆ ) exp ( (cid:219) M m ∆ ) + exp ( (cid:219) M ∆ ) , (3)where α L = . , α H = . are the lower and upper boundsof α and (cid:219) M m = . × − M (cid:12) yr − and ∆ = ( . / . ) × − M (cid:12) yr − are the midpoint and the width of the crossoverbetween α L and α H .For the case of hot accretion we have implemented amodule which calculates L inacc and distributes the energyacross the outer layers of the protostar. Each cell in thedeposition region receives the same specific energy additionsuch that the energy is distributed evenly with respect tothe stellar mass. For numerical stability, we set an upperlimit on the specific energy addition in each cell, such that ahigher L inacc results in distribution across a larger number ofcells and thus a deeper deposition of energy. We have testedthe robustness of this implementation by varying the energylimit for each cell so that the energy is distributed acrosslarger and smaller regions and found almost no variation inthe evolution of the protostars, given that convection in the surface layers in most cases readily redistribute the addedheat. In this section we present the stellar evolution models, theirtracks in the HR diagram, and compare models with differ-ent thermal efficiencies.
In Fig. 5 is shown the evolutionary tracks for a protostarwith a final mass M = . (cid:12) for six different thermal effi-ciencies. The accretion rate of the protostar is shown in thetop row of Fig. 4. The figure includes indicators for the igni-tion of deuterium as well as the position of the protostar at0.1 Myr and 1 Myr. We note that many of the stars remaindeeply embedded for the first ∼
500 kyr and therefore the ini-tial part of the tracks presented here should not be compareddirectly with the observations. Below we use a simple crite-ria to decide when a star is embedded based on a comparisonbetween the stellar and envelope mass based on the resultsof Frimann et al. (2016a). For a more refined estimate, a syn-thetic observation would have to be performed on the simu-lation, propagating the luminosity calculated from the stel-lar model out through the cloud environment (Frimann et al.2016a,b; Kuffmeier et al. 2017). If we disregard anisotropicshadowing and scattering of the light in the disk and en-velope the luminosity is (isotropically) preserved, while thespectral energy distribution and effective temperature is re-distributed to longer wavelengths. We have not added theaccretion luminosity to the bolometric luminosity in Fig.5 to highlight the differences in the protostellar evolution.Including L outacc reduces the relative differences between thetracks as the accretion luminosity dominates the protostellarluminosity during periods of high accretion rates.It is evident that the initial evolution of the protostaris highly sensitive to the thermal efficiency α . Models witha non-zero thermal efficiency above . undergo a rapid ini-tial expansion phase, where the hydrostatic structure of thestar is readjusted due to the large deposition of accretion en-ergy during the initial phase where accretion rates are high.This effect makes the evolution in hot models insensitive tothe exact size and mass choice of the initial second Larsoncore as the entropy in the star is largely dominated by theaccreted matter, and the hydrostatic equilibrium becomesindependent of the initial characteristics, only depending onthe thermal efficiency and the initial accretion rate. For mod-els with α < . the addition of matter leads to a contractionof the protostar and they quickly reach central temperatureshigh enough for the ignition of deuterium and subsequentlyhydrogen. The readjustment leads to a strong dependenceon the initial conditions for the cold models as discussedby Hosokawa et al. (2011). Cold models with a larger initialcore will not contract enough to reach the Zero-Age Main Se-quence (ZAMS) before accretion ceases and thus ignition ofhydrogen is delayed compared to the cold models presentedhere, which ignite hydrogen after roughly kyr while ac-cretion is still ongoing. The early ignition of hydrogen meansthat the models reach the ZAMS early on and subsequently MNRAS000
500 kyr and therefore the ini-tial part of the tracks presented here should not be compareddirectly with the observations. Below we use a simple crite-ria to decide when a star is embedded based on a comparisonbetween the stellar and envelope mass based on the resultsof Frimann et al. (2016a). For a more refined estimate, a syn-thetic observation would have to be performed on the simu-lation, propagating the luminosity calculated from the stel-lar model out through the cloud environment (Frimann et al.2016a,b; Kuffmeier et al. 2017). If we disregard anisotropicshadowing and scattering of the light in the disk and en-velope the luminosity is (isotropically) preserved, while thespectral energy distribution and effective temperature is re-distributed to longer wavelengths. We have not added theaccretion luminosity to the bolometric luminosity in Fig.5 to highlight the differences in the protostellar evolution.Including L outacc reduces the relative differences between thetracks as the accretion luminosity dominates the protostellarluminosity during periods of high accretion rates.It is evident that the initial evolution of the protostaris highly sensitive to the thermal efficiency α . Models witha non-zero thermal efficiency above . undergo a rapid ini-tial expansion phase, where the hydrostatic structure of thestar is readjusted due to the large deposition of accretion en-ergy during the initial phase where accretion rates are high.This effect makes the evolution in hot models insensitive tothe exact size and mass choice of the initial second Larsoncore as the entropy in the star is largely dominated by theaccreted matter, and the hydrostatic equilibrium becomesindependent of the initial characteristics, only depending onthe thermal efficiency and the initial accretion rate. For mod-els with α < . the addition of matter leads to a contractionof the protostar and they quickly reach central temperatureshigh enough for the ignition of deuterium and subsequentlyhydrogen. The readjustment leads to a strong dependenceon the initial conditions for the cold models as discussedby Hosokawa et al. (2011). Cold models with a larger initialcore will not contract enough to reach the Zero-Age Main Se-quence (ZAMS) before accretion ceases and thus ignition ofhydrogen is delayed compared to the cold models presentedhere, which ignite hydrogen after roughly kyr while ac-cretion is still ongoing. The early ignition of hydrogen meansthat the models reach the ZAMS early on and subsequently MNRAS000 , 1–20 (2017) he luminosity spread in young clusters Figure 4.
Thumbnails for 12 different stars (one per row, continued on next page) with a mass between 0.7 and 0.8 M (cid:12) at a clusterage of 2.3 Myr. The numbers in parenthesis indicated in the rightmost plot correspond to the numbers in Figs. 1, 6, and 7. Each columndensity plot is based on a (10,000 AU) cut-out from the model. From left to right the panels show the gas distribution at the birthof the star, after 100 kyr, 500 kyr, 1 Myr, and at a cluster age of 2.3 Myr. The magenta star at the center of each panel indicates thestar position, red circles are stars with masses M < . (cid:12) , brown triangles are stars with . (cid:12) < M < . (cid:12) , while orange squaresindicate stars with M > . (cid:12) . The right-most panel shows the accretion rate. The panels illustrate the effect of different environmentsand the variability of accretion histories in a single mass bracket. Notice how the presence of binary and multiple systems in the modelnaturally induces rapidly varying accretion rates for some of the stars. cross it as the continued accretion causes the models to havesmaller radii than ZAMS stars of equivalent mass. Thus thetracks of the cold models cross below the ZAMS as they con-tinue to accrete mass and climb toward higher luminositiesand higher effective temperatures. Once the mass exceeds M ∼ . (cid:12) or accretion stops the tracks rejoin the ZAMSagain, where they will settle at the locations correspond-ing to their final mass. An example of such an evolutionarytrack is seen in Fig. 5. Our model features smaller initialcores than previous studies (Baraffe et al. 2012, 2017; Ku- MNRAS , 1–20 (2017)
Jensen & Haugbølle
Figure 4 – continued from previous page nitomo et al. 2017) of accreting models and as a result thetracks for the cold models are different than previous studiesas the protostars remain more compact.The cold accretion models presented here follow almostidentical tracks in the HR diagram independent of the ac-cretion rate and final mass when the accretion luminosityis not included. This behaviour is similar to the models ofHosokawa et al. (2011) albeit the tracks are different dueto different initial conditions. Including the accretion lumi-nosity leads to a spread in luminosities due to differences inaccretion rates, but the effect is diminished as accretion ratesdrop below − M (cid:12) yr − . We are not convinced that the evo-lution of these purely cold models with α = . throughout the entire accretion phase are representative of real proto-stars. Physically, back-of-the-envelope calculations suggestthat for accretion rates (cid:219) M > − M (cid:12) yr − some fraction ofthe accretion luminosity will heat the protostar (see Ap-pendix B in Baraffe et al. (2012)). Similarly, a recent numer-ical study of the accretion shock onto accreting protoplanetsindicate that cold accretion is unlikely for canonical valuesfor the planetary accretion process (Marleau et al. 2017).Due to the compact structure, the purely cold models lacka contraction phase at the end of accretion. This leads to abimodal feature as the cold accreting PMS stars are eitherin a high luminosity state – due to the large accretion lu-minosities – or in a low luminosity state lying at or close to MNRAS000
Figure 4 – continued from previous page nitomo et al. 2017) of accreting models and as a result thetracks for the cold models are different than previous studiesas the protostars remain more compact.The cold accretion models presented here follow almostidentical tracks in the HR diagram independent of the ac-cretion rate and final mass when the accretion luminosityis not included. This behaviour is similar to the models ofHosokawa et al. (2011) albeit the tracks are different dueto different initial conditions. Including the accretion lumi-nosity leads to a spread in luminosities due to differences inaccretion rates, but the effect is diminished as accretion ratesdrop below − M (cid:12) yr − . We are not convinced that the evo-lution of these purely cold models with α = . throughout the entire accretion phase are representative of real proto-stars. Physically, back-of-the-envelope calculations suggestthat for accretion rates (cid:219) M > − M (cid:12) yr − some fraction ofthe accretion luminosity will heat the protostar (see Ap-pendix B in Baraffe et al. (2012)). Similarly, a recent numer-ical study of the accretion shock onto accreting protoplanetsindicate that cold accretion is unlikely for canonical valuesfor the planetary accretion process (Marleau et al. 2017).Due to the compact structure, the purely cold models lacka contraction phase at the end of accretion. This leads to abimodal feature as the cold accreting PMS stars are eitherin a high luminosity state – due to the large accretion lu-minosities – or in a low luminosity state lying at or close to MNRAS000 , 1–20 (2017) he luminosity spread in young clusters the ZAMS. Such a bimodal feature does not fit observationsof young clusters, which show signs of a longer contractionphase towards the ZAMS. We stress that we do not excludethat possibility that α = . for large parts of the accretionprocess where accretion rates are below ∼ − M (cid:12) yr − andit is possible that a small number of very low-mass objectsnever reach accretion rates of this order and evolve throughcold accretion.The choice of thermal efficiency has another profoundeffect on the models presented here. Purely cold models arefully convective until the ignition of hydrogen at which pointdeuterium burning continues in a shell and the region be-tween the burning zones is radiative. Outside the deuteriumburning shell the protostars remain convective until theysettle on the ZAMS at the end of accretion .Models with non-zero α exhibit only partial convectionduring the protostellar accretion phase. In a very small re-gion in the center there exist a radiatively dominated core,while high accretion rates can increase temperatures in theoutermost layers of the star. This drives a small temporaryradiative zone in the middle of the star and a dynamical in-stability, as discussed below. That convection is not stablein the outer layers is not surprising given the addition ofentropy during hot accretion which can alter the tempera-ture gradient and thus inhibit the convective motion. Thiseffect is reproduced in two dimensional models of accretingprotostars (Geroux et al. 2016). Despite the appearance oftemporary radiative zones in the outer layers, accreted deu-terium and lithium will reach the burning region and oncethe star reaches the PMS these regions disappear, exceptwhen the star is in a high accretion state, as seen in Fig. 9.The radiative nature of the innermost core during hot ac-cretion is more surprising. The same behaviour was foundby Stahler et al. (1980a) and it leads to a temperature in-version at the boundary between the radiative core and theconvective exterior region, which subsequently results in ashell ignition of deuterium. Stahler et al. (1980a) explainedthe occurrence of a temperature inversion between the con-vective region and the radiative core by the convectively sta-ble structure of their initial hydrostatic core. While in ourmodels the initial core is convective at the beginning of thesimulation, the addition of accretion energy in the hot ac-cretion models changes the temperature gradient drasticallyfrom the onset of accretion, and in particular in the earlyphase. This effect causes convection to halt near the bound-ary of the initial core and the region remains convectivelystable until the end of the main accretion phase. Once accre-tion ceases the protostar quickly reaches the same convectivelayout as main sequence stars of equal mass. Fig. 9 show star from Fig. 4 during the main accretion phase after the on-set of deuterium fusion, but just after a period with highaccretion rates. The temperature inversion is clearly visiblein the high density region, where the star is radiative.In Figs. 6 and 7 we present the evolution of the 12 pro-tostars from Fig. 4 in the HR diagram to illustrate howthe evolution of protostars in the same mass bracket is af-fected by differences in the time-dependent accretion. The Stars which reach central temperatures sufficient for CNO to bethe predominant fusion process become radiative in the envelopeas classical theory prescribe.
HR diagram does not include the accretion luminosity L outacc to highlight the protostellar variations. The luminosity ofthe protostars vary by almost 1 dex during the initial 1 Myrof the stellar life but converge to similar tracks as they reachthe contraction phase and descend on the Hayashi track.Figure 8 show the evolution of the stellar variables T eff , R and the luminosity from the initial accretion phase untilthe ZAMS for star ; the same protostar as in Fig. 5. Theinitial expansion phase is clearly seen in the early stageswhile the effect of late accretion bursts on the evolution isevident from several bursts occurring between kyr and Myr, which halt the contraction and cause large fluctuationsin the observable variables.In Figs. 5, 7, and 8 we see instabilities near the contrac-tion tracks towards the end of accretion. These fluctuationsin luminosity and surface temperature are driven by the de-position of accretion energy. To investigate this in more de-tail, in Fig. 9 is shown the T − ρ phase-space plot for star . The blue dot indicates the extent of the region whereexternal heat from the luminosity at the accretion shockis deposited. While it appears like a large region in phase-space, it only corresponds to 0.16% of the mass and 7.8%of the radius. The increased temperature in the outer layerscreates an unstable radiative region, where the adiabatic co-efficient decreases below γ ad < / leading to a dynamicalinstability. This can be seen in Fig. 9 where the slope of thecurve drops near log ρ = − . The region contracts releasingenergy, which causes the luminosity and effective tempera-ture to increase, while the star settles to a new stable state.After the region has contracted the adiabatic coefficient risesand convection is restored. A Similar instability is seen infigure 3 of Hosokawa et al. (2011) (case mE-C) for a modelfeaturing episodic accretion with a cold accretion shock andhigh accretion rates of (cid:219) M = − M (cid:12) yr − . As a consequence,in their models, a large amount of matter is added to theouter layers of the star with the same entropy as the existingphotosphere. The star responds by increasing the luminosityand effective temperature to radiate the additional energyand reestablish hydrostatic equilibrium. The protostar doesnot expand since the accretion timescale, t acc = M (cid:219) M , is shorterthan the thermal timescale, t KH = GM RL , which means thethermal equilibrium is not maintained. This is the case oncethe protostars reaches a certain mass for realistic accretionrates. The importance of the ratio between t KH and t acc foraccreting stars is also discussed by Prialnik & Livio (1985).To assess how robust the instability is to how energy is de-posited during hot accretion we have evolved models usingsmaller and larger regions for the energy deposition. Whilethis changes the details of where exactly the discontinuityoccurs and the early time evolution in a burst, it does notchange the overall behaviour. Even though both Hosokawaet al. (2011) and ourselves observe these instabilities, wecaution that the results are based on one-dimensional stel-lar evolution codes. A proper stability analysis is needed toconfirm and understand the nature of the instability. Gerouxet al. (2016) studied dynamic 2D models of accreting proto-stars and found that hot accretion changes and sometimeseven inverts the temperature gradient locally. This inhibitsconvection and leads to local dynamical instabilities as theones described above, however it is not clear to what extentthese effects are comparable. MNRAS , 1–20 (2017) Jensen & Haugbølle . . . . . . . eff [K] − . − . − . − . − . . . . l og L / L (cid:12) α = 0 . α = 0 . α = 0 . α = 0 . α = 0 . α ( ˙ M )ZAMS1M (cid:12) H ignition0.1 Myr1.0 Myr
Figure 5.
Models with varying thermal efficiencies for star inFig. 4 with M = . (cid:12) . The outward component of the accretionluminosity L outacc is not included. The environment and accretionprofile for this star is shown in the top row of Fig. 4. . . . . . . eff [K] − . − . − . − . . . . l og L / L (cid:12) ZAMS1M (cid:12) H ignition0.1 Myr1.0 Myr
Figure 6.
Evolutionary tracks for 5 stars with masses in the inter-val M ∈ [ . , . ] M (cid:12) with the dynamical thermal efficiency α ( (cid:219) M ) .The stars in this figure does not feature significant accretion burstin the PMS phase, opposed to the stars in Fig. 7. The outwardcomponent of the accretion luminosity L outacc is not included. Thenumbers refer to Fig. 4. The existence and strength of the instability vary signif-icantly between different protostars, depending on the inter-mittency of their accretionary history, as some hardly expe-rience the instability while it dominates the PMS evolutionfor others as seen in Figs. 5–7. The period of the instabili-ties is coupled to the change in accretion rates and can beanywhere from tens of yr to a kyr. They are therefore notsimilar to the observed pulsations in PMS stars, which areassumed to result from the κ and γ mechanisms or convec-tive barriers (see e.g. Zwintz 2015, and references therein).As discussed below, periodically changing accretion rates inbinary systems can lead to periodic changes in the total lu-minosity by more than a factor of 10 on a relatively shorttime-scale (see Fig. 14). . . . . . . eff [K] − . − . − . − . . . . l og L / L (cid:12) ZAMS1M (cid:12) H ignition0.1 Myr1.0 Myr
Figure 7.
Evolutionary tracks for 7 stars with masses in theinterval M ∈ [ . , . ] M (cid:12) with the dynamical thermal efficiency α ( (cid:219) M ) . The stars in this figure have significant accretion burststhat cause instabilities to occur. The outward component of theaccretion luminosity L outacc is not included. The numbers refer toFig. 4. Collinder 69 is a young open cluster located in the Orion con-stellation and part of the or λ -Orionis star forming region.It has been speculated (Cunha & Smith 1996) that a recentsupernova, less than 1 Myr ago, has exploded in Collinder69, quenching star formation in the cluster and creating amolecular shell at around 15 pc from the center of Collinder69 (Lang et al. 2000), where ongoing star formation stillcontinues today in cold clouds such as B30 and B35 (Koeniget al. 2015; Liu et al. 2016; Hu´elamo et al. 2017). The youngstellar objects (YSOs) in Collinder 69 includes class II andclass III stars, and stars with debris disks. The stellar pop-ulation in Collinder 69 is therefore more evolved than thestellar population in our model, which also contains new-born protostars. We have compensated for this by removingall embedded objects from our sample, as discussed below.It is estimated to be approximately 5 Myr old Bayo et al.(2011), though this estimate is based on classical isochrones,which could lead to significant errors, and Bell et al. (2013)find it to be up to 10 Myr old. Nonetheless, for simplicityand since this is the first paper to explore a self-consistentmodel of stellar evolution in an evolving molecular cloud, wewill consider this age to be the most likely for the cluster andcompare our models at 5 Myr with the inferred locations inthe HR diagram of the members of Collinder 69 from Bayoet al. (2011).The bolometric luminosities plotted in the HR diagramsare the total luminosities, which are the sum of the stellarluminosities and the outward component of the accretionluminosities L tot = L ∗ + L outacc . In the observational censusesall the embedded protostars will be absent, and we there-fore have to make an estimate for when a star is a Class0 protostar. In principle one should do a careful evaluationof the angle-averaged dust extinction for each star, or evencompute HR diagrams from many different viewing anglesof the cluster, and average the diagrams themselves. In thispaper, we do not try to make a detailed comparison between MNRAS000
Evolutionary tracks for 7 stars with masses in theinterval M ∈ [ . , . ] M (cid:12) with the dynamical thermal efficiency α ( (cid:219) M ) . The stars in this figure have significant accretion burststhat cause instabilities to occur. The outward component of theaccretion luminosity L outacc is not included. The numbers refer toFig. 4. Collinder 69 is a young open cluster located in the Orion con-stellation and part of the or λ -Orionis star forming region.It has been speculated (Cunha & Smith 1996) that a recentsupernova, less than 1 Myr ago, has exploded in Collinder69, quenching star formation in the cluster and creating amolecular shell at around 15 pc from the center of Collinder69 (Lang et al. 2000), where ongoing star formation stillcontinues today in cold clouds such as B30 and B35 (Koeniget al. 2015; Liu et al. 2016; Hu´elamo et al. 2017). The youngstellar objects (YSOs) in Collinder 69 includes class II andclass III stars, and stars with debris disks. The stellar pop-ulation in Collinder 69 is therefore more evolved than thestellar population in our model, which also contains new-born protostars. We have compensated for this by removingall embedded objects from our sample, as discussed below.It is estimated to be approximately 5 Myr old Bayo et al.(2011), though this estimate is based on classical isochrones,which could lead to significant errors, and Bell et al. (2013)find it to be up to 10 Myr old. Nonetheless, for simplicityand since this is the first paper to explore a self-consistentmodel of stellar evolution in an evolving molecular cloud, wewill consider this age to be the most likely for the cluster andcompare our models at 5 Myr with the inferred locations inthe HR diagram of the members of Collinder 69 from Bayoet al. (2011).The bolometric luminosities plotted in the HR diagramsare the total luminosities, which are the sum of the stellarluminosities and the outward component of the accretionluminosities L tot = L ∗ + L outacc . In the observational censusesall the embedded protostars will be absent, and we there-fore have to make an estimate for when a star is a Class0 protostar. In principle one should do a careful evaluationof the angle-averaged dust extinction for each star, or evencompute HR diagrams from many different viewing anglesof the cluster, and average the diagrams themselves. In thispaper, we do not try to make a detailed comparison between MNRAS000 , 1–20 (2017) he luminosity spread in young clusters Figure 8.
The time evolution of stellar variables T eff , R and the luminosity including the individual components L ∗ and L acc for star in Fig. 4 with the dynamical thermal efficiency α ( (cid:219) M ) . − − − ρ l og T (cid:15) F / kT ≈ rad ≈ P gas Age : .
74 kyrM : .
36 M (cid:12) M final : .
78 M (cid:12) ˙M : . × − M (cid:12) yr − T- ρ profile ConvectionNo mixingOvershoot
Burning limits:
H burning H burning Figure 9. T- ρ figure of star in Fig. 4. Colours indicate howheat is transported and matter is mixed. Fusion limits are indi-cated by dash dotted lines. The blue point indicate the depth ofthe region where accretion energy is distributed at that specificpoint in time. The star was evolved using the dynamical α ( (cid:219) M ) model. the modelled and observed clusters, and we have thereforedecided to use a simpler approach based on the correlationbetween envelope-stellar mass ratio and protostellar class inFrimann et al. (2016a). We define the envelope mass as thetotal gas mass inside a 10,000 AU distance from each starand compare it with the stellar mass. If the envelope massis larger than the stellar mass the star is assumed to embed- ded. Applying this selection criteria to the last snapshot inthe model, we find that 182 stars are embedded at 5 Myr,leaving 413 stars from the evolved population together with231 stars from the young population to be included in theHR diagrams.Fig. 10 show an estimate of the HR diagram for PMSstars with masses in the range M ∈ [ .
08; 1 . ] M (cid:12) in thesynthetic cluster at an absolute age of 5 Myr. We samplethe stellar distribution in a time interval of 200 kyr around5 Myr, and construct a 2D weighted histogram contour toshow the extent of the synthetic cluster in the HR diagram.We use a 200 kyr interval instead of a snapshot, to properlycapture, in a statistical sense, any episodic accretion eventsthat may change the luminosity significantly. The weightscorresponds to the time the star spend at that specific lo-cation in the HR diagram, and the histogram is therefore afair comparison with observations. The upper mass limit isbecause we use a mass range similar to the estimated massrange of the confirmed members of Collinder 69. We do notinclude brown dwarfs in the synthetic cluster even thoughthey are included in the Collinder 69 data as the mesa mod-els were not setup to handle the evolution of these, and theinitial mass function in the simulation may be incompleteat these masses. The figure includes Myr of star forma-tion stitched together from the . Myr molecular cloudsimulation as described in section 2.1.Fig. 11 is similar to Fig. 10 but features members of theOrion Nebular Cluster (ONC) and the mass is no longer lim-ited to . (cid:12) . For the ONC members there is a slight offsettowards higher luminosities between the high density regionof the synthetic cluster and the high density region of ONCmembers. This could be the result of several distinct stellarpopulations including a recent star formation burst (Beccari MNRAS , 1–20 (2017) Jensen & Haugbølle − − − l og L / L (cid:12) α = 0.0 (A) α = 0.01 (B) . . . log T eff [K] − − − l og L / L (cid:12) α = 0.25 (C) . . . log T eff [K] α ( ˙ M ) (D) ZAMSCollinder 691 Myr10 Myr − − − − − − − − − − − − − − − − − Figure 10.
2D histogram of the stellar positions in the HR diagram for the thermal efficiencies α = , . , . , α ( (cid:219) M ) in the syntheticcluster at a cluster age of 5 Myr. The density includes each position of the stellar models in the time interval t cluster ∈ [ . , . ] Myr toaccount for the possible variations in the 200 kyr interval due to episodic accretion events. The histogram is properly weighted accordingto the time interval each star spends in the different states. Only stellar models with a final mass in the range M ∈ [ . , . ] M (cid:12) areincluded. Note that the histogram density is log-scaled. We caution that at first sight it is easy to overestimate e.g. the number of starsin the high-luminosity state. We have chosen this scaling to show the full extent of the distribution. Black points are confirmed membersof Collinder 69 (Bayo et al. 2011). 1 Myr and 10 Myr isochrones are non-accreting pre-main sequence models evolved in mesa with thesame composition. et al. 2017). This is not entirely reproduced in the syntheticcluster simulations, due to the limited simulation time spanand possibly also the extent of the simulated molecular cloudpatch. The IMF of the ONC also show a peaked distributionsomewhat different from the Chabrier IMF and the IMF inour model (Dib 2014). Younger stars created within the lastMyr would appear with larger luminosities in the syntheticcluster and could explain the discrepancy between the pop-ulations.Regarding Collinder 69 we note that the estimated ageof 5 Myr is based on the best fit to the observed dataamongst the cluster members with low T eff (Bayo et al.2011). The choice to estimate the age based on the coolercluster members is due to a smaller luminosity spread inthis part of the HR diagram. In the cooler part of the HRdiagram 95% of the cluster members lie above a 20 Myrisochrone while a larger age spread exist for the warmermembers. The mass range for the members of Collinder 69is inferred from classical isochrones assuming a cluster age of 5 Myr. Bayo et al. (2012) found the mass estimates to berobust under variations in the cluster age. We do not be-lieve this conclusion would change significantly in a frame-work of accreting PMS models due to similarities betweenthe Hayashi tracks of accreting models and non-accretingmodels of similar mass. The tracks of low- and high-massstars do coincide at early stages for accreting models but atthis stage the stars are expected to be embedded and shouldthus not distort the mass estimates. In a realistic model of a molecular cloud complex dominatedby turbulent fragmentation massive stars are not formedimmediately from as single large core, but instead is fedthrough the filaments and take a significant amount of timeto accrete, reaching several Myr for the most massive starsin the simulation (Padoan et al. 2014b). This changes ourpicture of stellar evolution. To illustrate it, in Fig. 13 are
MNRAS000
MNRAS000 , 1–20 (2017) he luminosity spread in young clusters examples of low-, intermediate- and high-mass stellar tracksusing the variable α ( (cid:219) M ) -model for the thermal efficiency to-gether with arrows indicating different phases of this moredynamic picture of early stellar evolution.Initially, for the first ∼
100 kyr all stars evolve in a sim-ilar manner, starting from a universal second Larson core(A). Following almost identical tracks in the HR diagramthey rapidly expand and increase their effective temperature(B). Towards the end of the main accretion period duringthe deeply embedded stage, low-mass stars start to contract,descending along the Hayashi track (C). If the star insteadkeeps accreting, the increase in temperature continues andwhen the core becomes radiative, and deuterium burningcontinues in a shell, it starts to expand rapidly. Eventually,the outer layers become radiative and the high-mass star be-gin to follow a Henyey track at a slightly inclined luminosity,due to the ongoing mass accretion. The (cid:12) star shown inFig. 13 reaches the Henyey track (D) at an age of ∼
800 kyrand a mass of ∼ . (cid:12) , and stays on it during ∼
300 kyr.Notice that the evolution of massive stars as they join theHenyey track is almost independent of the final mass of thestar, and more dependent on how much mass they are ableto accrete before the shell-burning of deuterium commence.In this case, once the mass reaches ∼ (cid:12) the star ignitesthe p-p chain in the core and soon after joins the ZAMS. Atthis stage the star is massive enough that the thermal effi-ciency of the accretion has insignificant effect on the stellarstructure and evolution. Finally, the high-mass star in Fig.13 during a period of ∼
500 kyr moves on the ZAMS, almostdoubling the mass and reaching a final mass of 8 M (cid:12) (E).It is an important conceptual change compared to theclassic picture of non-accreting stellar evolution segmentedin to low-mass Hayashi tracks and high-mass Henyey tracks.While some details may change when the effects of e.g. rota-tion, magnetic fields, outflows, and radiation are eventuallyincluded in the models, the evolution is fundamentally aconsequence of the time-scales necessary for accreting themass, and therefore robust. The corresponding upper enve-lope can be appreciated in Fig. 11. While our tracks superfi-cially are similar to those of Haemmerl´e et al. (2013, 2016),the difference is the time-scales, since they consider rapidcore-accretion, which is not supported by our global modelsof star forming clouds.Exactly at which mass the high-mass stars start to movealong the Henyey track depends on the early accretion rate,which is regulated by the rate of infalling material fromlarger scales. In the case of a free-fall isothermal envelopethis is bounded by ∼ c s / G , where the exact numerical pre-factor is between 1 and 50 (Penston 1969; Larson 1969; Shu1977; Hunter 1977; Foster & Chevalier 1993), and thereforedepends on the temperature of the envelope and surroundingmolecular cloud core.The above description of the evolution applies to hot ac-cretion models. As already discussed, cold accretion modelsbehave very differently. Because of the low initial tempera-ture, the stars rapidly contract, heat up, and ignite hydro-gen. Even low-mass stars reach the ZAMS as early as after100 kyr, and then start to evolve along the ZAMS whileaccreting. The late evolution, in hot accretion models, forhigh-mass stars along the ZAMS is basically for the samereasons: at the point where the star reaches the ZAMS thetotal luminosity is completely dominated by radiation from the surface, while the accretion luminosity is sub-dominant,and the evolution proceeds as if it was cold accretion. Interactions between stars on small scales and between starsand the collective gas and stellar potential on larger scalescan affect the accretion history and thereby the stellar evo-lution.On the larger molecular cloud scale, as protostars movearound in the collective gravitational potential, they mayenter or exit high density regions and filaments, which cancause the time evolution in the accretion rates to deviatedrastically from models of isolated cores. In Figs. 1 and 4we see that some stars cluster together in dense regions whileothers evolve quietly in isolation (see also Kuffmeier et al.2017). The accretion profiles show large variations with ageneral trend where unsteady accretion rates are associatedwith stars in dense regions harboured in multiple systems,while isolated stars have smoother and monotonically de-creasing accretion profiles.Approximately half of all young solar mass stars arefound in binary or multiple systems (Tobin et al. 2016).This can induce periodically varying accretion rates (Padoanet al. 2014b) and affect the protostellar and PMS evolution.In Figs 14 and 15 we show an example of a protostar with abinary induced periodicity in the accretion rate. At an ageof 900 kyr the shown star has a mass of 0.5 M (cid:12) , while itsprimary companion is a 1.3 M (cid:12) star. The system has a semi-major axis of a = . AU and an eccentricity of e = . .The period is 63.5 yr. This periodicity gives large fluctua-tions in the effective temperature and halts the contractionof the protostar, thus delaying the evolution towards theZAMS. The periodic accretion bursts also lead to luminos-ity fluctuations of more than one order of magnitude in syncwith the orbit time. This is because the secondary is on anelliptic orbit, where at aphelion it has the smallest veloc-ity difference with the envelope, and therefore the largestaccretion rate.As a result it is not readily evident to an observer atwhat evolutionary stage the star is in solely based on the ac-cretion rate, and a binary system can have luminosity varia-tions on the same scale as FU Orionis-like bursts with time-scales of tens to hundreds of years and periodic in nature.Such bursting binary systems could occur both in the proto-stellar and the PMS evolutionary phases. Recently, severaldeeply embedded binary protostars have been found with arelatively small seperation and clear chemical fingerprints ofrecent episodic accretion events (Frimann et al. 2017). One motivation behind the development of accreting proto-stellar models have been to address the observed luminosityspread in young stellar clusters as proposed by Baraffe et al.(2009). Figs 10 and 11 show good agreement between the lu-minosity spread in the observed HR diagrams of two young
MNRAS , 1–20 (2017) Jensen & Haugbølle − − l og L / L (cid:12) α = 0.0 (A) α = 0.01 (B) . . . . . log T eff [K] − − l og L / L (cid:12) α = 0.25 (C) . . . . . log T eff [K] α ( ˙ M ) (D) ZAMSONC − − − − − − − − − − − − − − − − − Figure 11.
2D histogram of the stellar positions in the HR diagram for the thermal efficiencies α = , . , . , α ( (cid:219) M ) in the syntheticcluster at a cluster age of 5 Myr. The density includes each position of the stellar models in the time interval t cluster ∈ [ . , . ] Myr toaccount for the possible variations in the 200 kyr interval due to episodic accretion events. The positions are time-weighted. Note thatthe histogram density is log-scaled. We caution that at first sight it is easy to overestimate e.g. the number of stars in the high-luminositystate. We have chosen this scaling to show the full extent of the distribution. Black points are confirmed members of the Orion NebularCluster (ONC) (Da Rio et al. 2010). Density is given in counts per bin. clusters and the luminosity spread in the synthetic cluster ifone allows some variation in α values in a coeval population.The luminosity spread in the synthetic cluster is not solelya consequence of the episodic nature of the accretion pro-files or the implementation of the thermal efficiency of theaccretion shock, but is the result of a number of differentfactors.The variations in the thermal efficiency of the accretionshock produce large variations in the protostellar evolution-ary tracks, but the effect is mostly limited to the proto-stellar phase ( t star < Myr), where variations in accretionrates in general are most pronounced causing large fluctua-tions in the luminosity and radii of the individual models.Once the protostars enter the PMS characterised by low ac-cretion rates with occasional accretion bursts, the variationsare much less pronounced. For the warm models with α > . the tracks for different thermal efficiencies converge towardseach other and the memory of the initial accretion historyand the differences in thermal efficiency is essentially erased.From Fig. 5 we see that these models reach the Myr pointat almost identical positions in the HR diagram and follow similar tracks towards the ZAMS. The PMS tracks of thewarm models are very similar to the classical Hayashi tracksstars with the same masses. In the purely cold models with α = . stars follow very different tracks and have much lowerluminosities, however as discussed in section 3.1 we do notbelieve the evolution of the purely cold models, with α = . throughout the entire accretion phase, are representative ofthe majority of real protostars. In Fig. 5 we see that themodel with α = . has a much shorter contraction phasethan the warmer models. Models with a thermal efficiencyin the range . < α < . can produce a substantial spreadin luminosity as the length of the contraction tracks and thetemporal evolution differ. Thermal efficiencies in this rangeare necessary to explain the faintest members in Collinder69 as seen in panels A and B in Fig. 10. The shorter contrac-tion phase of the stellar models with low thermal efficiencymeans that these would appear older in the old paradigm ofnon-accreting isochrones, where the relatively low luminosityof these models would not be reached until late in the PMSphase. Comparing the non-accretion Myr isochrone fromBaraffe et al. (2015) with the lower boundary of the accret-
MNRAS000
MNRAS000 , 1–20 (2017) he luminosity spread in young clusters − − − l og L / L (cid:12) classical (A) ˙ M = 10 − M (cid:12) yr − α = 0.25 (B) . . . log T eff [K] − − − l og L / L (cid:12) ˙ M = 10 − M (cid:12) yr − α = 0.01 (C) . . . log T eff [K] α ( ˙ M )no age spread (D) ZAMSCollinder 691 Myr10 Myr − − − − − − − − − − − − − − − − Figure 12.
Similar to Fig. 10 but with different models. Panel A: classical PMS models evolved in mesa but with the same formationhistory as the synthetic cluster models. Panel B: constant accretion rate (cid:219) M = − M (cid:12) yr − and α = . , again with the same formationhistory as the synthetic cluster. Panel C: constant accretion rate (cid:219) M = − M (cid:12) yr − and α = . , again with the same formation historyas the synthetic cluster. Panel D: the same as panel D in Fig. 10 but without the age spread of the formation history from the syntheticcluster, meaning that all stars have evolved for 5 Myr and stopped accreting. The high luminosity tracks in panel B and C are modelswhich are still actively accreting. ing models with non-zero α in Fig. 10 the accreting modelswould appear much older if the non-accreting isochrones wasused to estimate the age.Allowing for significant variations in thermal efficienciesbetween coeval stellar populations can produce an additionalluminosity spread. In figure 10 a combination of panels B andD in would almost reproduce the entire Collinder 69 lumi-nosity spread. However this combination requires that somestars have a low thermal efficiency throughout the entire evo-lution including the initial main accretion phase where accre-tion rates usually exceed (cid:219) M = − M (cid:12) yr − . It is uncertainif this is physically feasable. The purpose of the dynamical α model was to replicate a combination of several α values,however our model does not suffice in this regard as all theprotostars undergo hot accretion at some stage. It should benoted that the maximum thermal efficiency presented here, (cid:15) = . & α = . , is a rather conservative upper limit for thethermal efficiency compared with the theoretical maximumefficiency of from Stahler et al. (1980a). Additionally,we do consider it likely that some of the observed clustermembers in figs 10 and 11 are unresolved binaries, which will have a higher apparent luminosity than the individualstars.The introduction of realistic accretion histories createlarge variations in the evolutionary tracks in the early stagesof accretion, but the effects are again limited after t star > Myr. It should be noted though, that due to the limitedresolution in our model (50 AU), we are only capturing theinfall rate of material from the envelope to the disk, andusing it as a proxy for the the accretion rate of material fromthe disk to the star. Given sufficient numerical resolution andrelevant microphysics in the model to resolve and correctlyevolve the disks around the stars, a higher intermittency intime and large variability in amplitude of the accretion ratesare to be expected (Kuffmeier et al. 2017).To construct a model that captures the probable differ-ences in the accretion flow as a function of the accretion rate,we have introduced the dynamic thermal efficiency α ( (cid:219) M ) ,which further enhances the significance of the accretion his-tories. In Fig. 6 and Fig. 7 we illustrate differences in theevolution tracks of stars with the dynamical thermal effi-ciency and varying accretion rates. We see that there exists MNRAS , 1–20 (2017) Jensen & Haugbølle
Figure 13.
HR diagram with evolutionary tracks for stars of , , and (cid:12) with α ( (cid:219) M ) . The dash dotted line (dark blue, pink,and purple) show the protostellar luminosity L ∗ while the fullline (blue, orange, and green) shows the total luminosity L total = L ∗ + L acc . (A) is the initial expansion phase as the accretionbegins from the universal initial condition. (B) marks the mainaccretion phase, which is independent of the final mass of the star.(C) marks the Hayashi contraction track for low-mass stars. (D)marks the Henyey tracks for intermediate- and high-mass starswith continued accretion during hydrogen fusion. (E) marks thefinal accretion phase for high-mass stars along the ZAMS. T e ff [ K ] − L u m i n o s i t y [ L (cid:12) ] L tot L acc L star time since star formation [yr] − − − ˙ M [ M (cid:12) y r − ] R a d i u s [ R (cid:12) ] Figure 14.
Example of a protostar with periodic accretion ratesinduced by binary interactions. Top panels shows the evolution of T eff and R . Middle panel show the total luminosity as well as thetwo component L acc and L star . Lower panel shows the accretionrate. Fig. 15 show a smaller timespan where the periodicity isclearly visible. T e ff [ K ] − L u m i n o s i t y [ L (cid:12) ] L tot L acc L star − − − ˙ M [ M (cid:12) y r − ] . . . . . R a d i u s [ R (cid:12) ] Figure 15.
Same panel as Fig. 14 but with a closeup over a 1000yr period highlighting the periodicity in the accretion rate. minor temporal variations between the models as they reachthe PMS contraction tracks, with a ∼ . dex spread of the Myr markers. These variations go some way towards ex-plaining the observed luminosity spread in young clusters,but overall the PMS evolution of the models with the samethermal efficiency is not that different and the effect of vary-ing accretion histories is not sufficient to account for thewhole luminosity spread.From the global molecular cloud simulations we obtaina realistic age spread of protostars in a coeval population,something which has not previously been included in stud-ies of accreting protostellar models in relation to the stellarstructure and reproducing the luminosity spread in the HRdiagram. Including the absolute formation time of each pro-tostar naturally increases the spread in luminosities as somestars are much younger than the age of the cluster itself.This effect is important for replicating the observed lumi-nosity spread.In Fig. 12 panels A,B, and C we show classical PMSmodels and constant accretion models evolved in mesa fol-lowing the same birth times and stellar masses as the pro-tostars of the global molecular cloud models. Panel D showsthe dynamical α model with the realistic accretion rates butexcluding the intrinsic age spread from the cluster simula-tion. Removing the age spread all the stars are presented atan age fo ± . Myr. At this stage accretion is ceased and,as discussed above, the luminosity variations caused by thedifferent accretion profiles is significantly diminished. Thisillustrates the importance of the natural age spread of thecoeval stellar populations. When the natural age spread isconsidered some stars are actively accreting at different rates
MNRAS000
MNRAS000 , 1–20 (2017) he luminosity spread in young clusters while others are past the accretion phase and this effect nat-urally contributes to a significant luminosity spread. Boththe classical and the constant accretion models produce adecent luminosity spread when the stellar age spread is in-cluded. Comparing figs. 10 and 12 we see that ultimately therealistic accretion rates does increase the luminosity spreadand fit the observed data better. Specifically comparing thepanel B in Fig. 10 with panel C in Fig. 12 the realistic accre-tion rates produce more faint objects whereas that modelswith similar α but constant accretion rate fail to produceobjects with luminosities similar to the faintest objects inCollinder 69. The variable accretion rates also produce morehigh luminosity objects at a cloud age of Myr which is dueto the longer accretion timescales of the accretion profilesobtained from the molecular cloud simulation.The global models also introduce intrinsic effects of theenvironment to the accretion rates. The local density andgravitational perturbations between protostars in star form-ing regions provide variations in accretion histories, whichcan provide fluctuations in luminosity beyond the main ac-cretion phase. Such variations are evident in figure Fig. 4,where the extend of the accretion profiles vary significantlybetween different protostars of similar final mass. Some ex-hibit large accretion burst beyond the Myr mark, whileothers have finished accretion in less than ∼ kyr. Binaryinteractions are also visible for some of the accretion profilesas discussed in section 3.4. Accounting for these variationsrequires global simulations with high dynamical range inboth space and time. In particular, secondary stars in bi-nary systems on elliptic orbits will preferentially accrete thegas when close to the aphelion, where they spend most ofthe time and the velocity differential with respect to the gasis the smallest. This induces regular periodic and sometimesdramatic changes in the accretion rate (see also Padoan et al.2014b). In the case of short-period orbits in a deeply em-bedded core where accretion rates can reach ∼ − M (cid:12) yr − in the actively accreting state, this “forcing” could inducevariability in the stellar structure with periods equal to thebinary period creating an externally forced variable star asseen in figs. 14 and 15.Kuffmeier et al. (2017) presents zoom-simulations for 6protostars from a similar molecular cloud simulation with anincreased spatial resolution down to 0.6 AU. At this resolu-tion the accretion disk and accretion bursts driven by diskinstabilities are better resolved. These accretion burst canfurther increase the luminosity spread of the accreting pro-tostellar models and increase the dependence on the thermalefficiency.In recent works, Baraffe et al. (2017) and Vorobyovet al. (2017) have studied accreting PMS stars with accretionhistories from two-dimensional hydrodynamical disk simula-tions, which feature accretion burst driven by gravitationalinstabilities within the disk. These models complement ourmodels as they are able to resolve disk instabilities betterthan the global molecular cloud simulations utilised here.However they lack the additional effects from the local cloudenvironment, the ability to track the absolute cluster ageand both ideal and non-ideal magnetohydrodynamical ef-fects, which could alter the circumstellar disk structure andthe characteristics of the accretion bursts. Both groups alsoexperiment with varying thermal efficiencies as well as hy-brid α models similar to the α ( (cid:219) M ) model presented here. Baraffe et al. (2017) find a bimodal luminosity distributionin their HR diagrams with the cold models producing veryfaint objects and the hybrid models producing brighter ob-jects which lie close to the classical isochrones. To repro-duce the luminosity spread they suggest it is necessary toallow varying values of α . These conclusions are in accor-dance with our findings where a combination of panel Band D in Fig. 10 would produce the best fit to Collinder69. Kunitomo et al. (2017) recently investigated the effectsand implications of varying deuterium mass fractions in ac-creting models of PMS stars. They found that variationsin the deuterium mass fraction have insignificant effects onhot accreting models, while cold or slightly warm accretionmodels can exhibit large evolutionary variations dependingon the deuterium abundance. To which degree larger varia-tions in abundances within a cluster is feasible depends onthe amount of mixing and pollution by stellar winds andexploding stars in the environment (Vasileiadis et al. 2013).In summary, the HR diagram of the synthetic clusterpopulations are in good agreement with the HR diagrams ofobserved young cluster as is seen in Figs. 10 and 11, espe-cially if we allow some protostars to accrete with a low ther-mal efficiency throughout the entire formation. This agree-ment would most probably only improve, if we were ableto capture the additionally variability in the accretion rateinduced by circumstellar disks in the model.Developing isochrones to improve the age estimation ofyoung clusters and PMS stars is beyond the scope of thispaper, but we suggest the following method as a viable wayto improve cluster age determinations in future work. First,a stellar population from the synthetic cluster simulationis sampled N times for a specific cluster age. The stellarpopulation should be similar to that of the observationaldata, e.g the mass should be in the correct mass bracketand number of stars should be the same. Subsequently, eachsample is compared to the observed positions in the HR di-agram using a Kolmogorov-Smirnov (KS) test to determineto what confidence level we can reject the null-hypothesisthat the underlying distributions match. Choosing a reason-able confidence level, the number of rejections is calculatedand used to estimate the whether the synthetic cluster ageis a good match for the observed cluster population. Thisprocedure is repeated for a number of synthetic cluster agesand the cluster age with the least rejections is consider tobe the most likely. Such a procedure could provide a solidage estimates for young clusters, but the choice of thermalefficiency would remain an free parameter, something whichfuture work would need to address. A similar statistical ap-proach could also be developed to estimate the mass of PMSstars in young clusters, and the probability that they are in ahigh-luminosity state. This could be done using the distribu-tion of stellar models that are contributing to model pixelsoverlapping with the observed star in the HR diagram. The development of realistic co-evolved models of stellar evo-lution is important not only in the context of star formationitself, but also in relation to the study of protoplanetary where N is chosen to be sufficiently largeMNRAS , 1–20 (2017) Jensen & Haugbølle disk structure, planet formation and the early solar systemas improved models of irradiation play a vital role in thephysics and chemistry of circumstellar disks. In Fig. 5 wesee that large variations in thermal efficiency α have a pro-found effect on the stellar luminosity as well as the effectiveradius and temperature, primarily in the early stages of starformation. Recent observations have shown that circumstel-lar disks appear while the protostar is still embedded in-creasing the importance of a realistic protostellar structuremodel for a proper account of photo chemistry, dust evolu-tion, and the general development of models of circumstellardisks (Murillo et al. 2013; Yen et al. 2017).We have compared observations of two young clusterswith a synthetic cluster with varying thermal efficiencies andfound that modest thermal efficiencies are sufficient to pro-duce comparable luminosity spreads in the synthetic cluster.Obtaining bounds on α with the models and data presentedhere remains difficult as uncertainties in observational meth-ods as well as simplifications in the numerical models causefine tuning of the efficiency to be inconclusive at this point.Looking at panel B in Fig. 10 the models with α = . fit theobserved data with low luminosities very well. Due to accre-tion bursts, higher luminosities can also be explained at thislow thermal efficiency, though the density of observed starsat higher luminosity suggest that some stars must have atleast some periods of higher thermal efficiency or be youngerthan the modelled population to reproduce the observed lu-minosity spread.The dynamical model α ( (cid:219) M ) gives a reasonable fit tothe observations and can be used as a phenomenological ap-proximation until better models for the thermal efficiencyare developed. It should be stressed that even though webelieve the model could be a good approximation given thegeneral understanding we have of protostellar accretion, itdoes not have a strong quantitative theoretical basis.As noted by Baraffe et al. (2017), constraining the ther-mal efficiency is a formidable problem of radiative hydro-dynamics, which neither theoretical nor observational ap-proaches have been able to solve yet. In a recent paper,Marleau et al. (2017) studied the problem in the context ofplanet formation using one-dimensional radiative hydrody-namics to model the shock above a hydrostatic planetary at-mosphere. Their model does not include the effects of higherdimensionality or magnetic fields, but provides a first steptowards solving the problem with numerical models. Assum-ing the dynamics is comparable to the protostellar case wecan use these results as a guideline for reasonable thermalefficiencies. They found that the thermal efficiency varies de-pending on accretion rates, planetary mass and shock radii,which shows that large variations in efficiency can occur de-pending on the specific properties of the protostar and thelocal environment. They also find that even if all the kineticenergy is effectively converted to radiation at the shock, theaccreting matter can trap a large fraction of the energy,which will then be reaccreted. If this is the case, purelycold accretion does not seem possible during periods withhigh accretion rates, which is in line with the conclusions wehave made based on our cold models. These results seems tofavour a model where high accretion rates result in higherthermal efficiencies as is the case with the dynamical modelpresented here. In this exploratory paper we have presented a new method-ology for determining ages of stellar associations by usingsynthetic clusters simulations combining accretion historiesfrom global molecular cloud simulations with accreting pro-tostellar structure models with varying thermal efficiency.The molecular cloud simulation has provided us with re-alistic accretion rates accounting for the dynamical natureof the star forming region as well as the complete star for-mation history within the cloud enabling a more realisticcomparison to observations of young clusters. Ongoing starformation in the molecular cloud simulation, in accordancewith what is seen observationally for the Lambda-Orionisstar forming region and the ONC, is important for a suffi-cient age spread in the stellar population, and the presenceof young stars that are still undergoing periods of increasedaccretion.We have experimented with a number of different ther-mal efficiencies including a dynamical model dependent onthe instantaneous accretion rate. We find that protostarsevolved with cold accretion throughout the entire accretionprocess are less likely to be compatible with observationsof young stellar objects, however a small fraction protostarsevolved with insignificant accretion energy are necessary toexplain the faintest observations. Meanwhile a range of non-zero thermal efficiencies are in good agreement with majorityof the observed cluster members. Hot accretion with thermalefficiencies above 0.5 does not appear to be required to ex-plain the observed HR diagrams of young clusters. Given themany competing factors that enters into creating a specificsynthetic HR diagram, constraining the thermal efficiencyfurther is not currently possible.We construct a dynamical model for the thermal effi-ciency depending on the accretion rate α ( (cid:219) M ) . This modelscannot explain the entire luminosity spread by itself, butin combination with models with low thermal efficiencies of α ∼ . the observed luminosity spread of the young clusterCollinder 69 is reproduced within the estimated Myr age.We also find reasonable agreement between the luminosityspread of the Orion Nebular Cluster and the synthetic clus-ter though the former show signs of several star bursts over alarger volume, which is not the case for the synthetic cluster.Extending the run time, size and mass of the global molec-ular cloud simulation model to be closer to that of the ONCcould possibly improve the concordance between syntheticmodels and the ONC.In binary systems where the secondary is on an ellipticorbit the infall from the envelope can induce variability inthe stellar structure with periods equal to the binary periodcreating an externally forced variable star. The changes intotal luminosity can be of up to two orders of magnitude,which is close to what is seen for both embedded protostarsand FU Orionis like stars.The large time-scales involved in massive star formationand the existence of a universal initial condition changes theevolution history of massive stars, with a relatively largetime-span first evolving diagonally in the HR diagram, thenon a slightly inclined Henyey-track due to the ongoing ac-cretion, and finally moving along the ZAMS, while reachingthe final mass.In conclusion, to obtain precise ages for young stellar
MNRAS000
MNRAS000 , 1–20 (2017) he luminosity spread in young clusters clusters we need to go beyond simple isochrones and con-sider the combined impact of (i) realistic, time-dependent,accretion rates, (ii) warm accretion and (iii) the formationtime of individual stars This can be done using stellar evo-lution models obtained from simulations of star forming re-gions. Only then the observed luminosity spread can be fullyaccounted for. The numerical data sets used in this publica-tion, including the stellar structure evolution and gas prop-erties in the vicinity of the protostars, are available uponrequest to the authors. ACKNOWLEDGEMENTS
We would like to thank Remo Collet, Søren Frimann, G¨unterHoudek, Jes Jørgensen, Michael Kuffmeier, Paolo Padoan,and Mads Sørensen for valuable comments and discussionswhile carrying out this work. We are grateful to the anony-mous referee for constructive comments that helped to im-prove the manuscript. TH was supported by a Sapere AudeStarting Grant from the Danish Council for IndependentResearch. Research at Centre for Star and Planet Forma-tion is funded by the Danish National Research Foundation(DNRF97). We acknowledge PRACE for awarding us accessto the computing resource CURIE based in France at CEAfor carrying out part of the simulations. Archival storage andcomputing nodes at the University of Copenhagen HPC cen-ter, funded with a research grant (VKR023406) from VillumFonden, were used for carrying out part of the simulationsand the post-processing.
REFERENCES
Baraffe I., Chabrier G., Gallardo J., 2009, ApJ, 702, L27Baraffe I., Vorobyov E., Chabrier G., 2012, ApJ, 756, 118Baraffe I., Homeier D., Allard F., Chabrier G., 2015, A&A, 577,A42Baraffe I., Elbakyan V. G., Vorobyov E. I., Chabrier G., 2017,A&A, 597, A19Bayo A., et al., 2011, A&A, 536, A63Bayo A., Barrado D., Hu´elamo N., Morales-Calder´on M., MeloC., Stauffer J., Stelzer B., 2012, A&A, 547, A80Beccari G., et al., 2017, preprint, ( arXiv:1705.09496 )Bell C. P. M., Naylor T., Mayne N. J., Jeffries R. D., LittlefairS. P., 2013, MNRAS, 434, 806Carr J. S., 2007, in Bouvier J., Appenzeller I., eds, IAU Sym-posium Vol. 243, Star-Disk Interaction in Young Stars. pp135–146, doi:10.1017/S1743921307009490Chabrier G., 2003, PASP, 115, 763Choi J., Dotter A., Conroy C., Cantiello M., Paxton B., JohnsonB. D., 2016, ApJ, 823, 102Cox J. P., Giuli R. T., 1968, Principles of stellar structureCrutcher R. M., 2012, ARA&A, 50, 29Cunha K., Smith V. V., 1996, A&A, 309, 892Da Rio N., Robberto M., Soderblom D. R., Panagia N., Hillen-brand L. A., Palla F., Stassun K. G., 2010, ApJ, 722, 1092Dib S., 2014, MNRAS, 444, 1957Dobbs C. L., et al., 2014, Protostars and Planets VI, pp 3–26Eisner J. A., et al., 2010, ApJ, 718, 774Evans II N. J., et al., 2009, ApJS, 181, 321Ferri`ere K. M., 2001, Reviews of Modern Physics, 73, 1031Foster P. N., Chevalier R. A., 1993, ApJ, 416, 303Freedman R. S., Marley M. S., Lodders K., 2008, ApJS, 174, 504Frimann S., Jørgensen J. K., Haugbølle T., 2016a, A&A, 587, A59 Frimann S., Jørgensen J. K., Padoan P., Haugbølle T., 2016b,A&A, 587, A60Frimann S., et al., 2017, preprint, ( arXiv:1703.10225 )Geroux C., et al., 2016, A&A, 588, A85Grevesse N., Sauval A. J., 1998, Space Sci. Rev., 85, 161Haemmerl´e L., Eggenberger P., Meynet G., Maeder A., Charbon-nel C., 2013, A&A, 557, A112Haemmerl´e L., Eggenberger P., Meynet G., Maeder A., Charbon-nel C., 2016, A&A, 585, A65Hartmann L., 2008, Accretion Processes in Star Formation, 2edn. Cambridge Astrophysics, Cambridge University Press,doi:10.1017/CBO9780511552090Hartmann L., Cassen P., Kenyon S. J., 1997, ApJ, 475, 770Hartmann L., Herczeg G., Calvet N., 2016, ARA&A, 54, 135Haugbølle T., Padoan P., Nordlund A., 2017, preprint,( arXiv:1709.01078 )Herbig G. H., 1966, Vistas in Astron., 8, 109Herbig G. H., 1977, ApJ, 217, 693Herbig G. H., 1989, ESO Workshop on Low Mass Star Formationand Pre-Main Sequence Objects, 33, 233Heyer M., Krawczyk C., Duval J., Jackson J. M., 2009, ApJ, 699,1092Hosokawa T., Offner S. S. R., Krumholz M. R., 2011, ApJ, 738,140Hu´elamo N., et al., 2017, A&A, 597, A17Hunter C., 1977, ApJ, 218, 834Hunter T. R., et al., 2017, ApJ, 837, L29Jørgensen J. K., Visser R., Williams J. P., Bergin E. A., 2015,A&A, 579, A23Kauffmann J., Pillai T., Goldsmith P. F., 2013, ApJ, 779, 185Kenyon S. J., Hartmann L. W., Strom K. M., Strom S. E., 1990,AJ, 99, 869Koenig X., Hillenbrand L. A., Padgett D. L., DeFelippis D., 2015,AJ, 150, 100K´osp´al ´A., ´Abrah´am P., Prusti T., Acosta-Pulido J., Hony S.,Mo´or A., Siebenmorgen R., 2007, A&A, 470, 211Kuffmeier M., Frimann S., Jensen S. S., Haugbølle T., 2017, sub-mitted to MNRASKunitomo M., Guillot T., Takeuchi T., Ida S., 2017, A&A, 599,A49Lang W. J., Masheder M. R. W., Dame T. M., Thaddeus P., 2000,A&A, 357, 1001Larson R. B., 1969, MNRAS, 145, 271Liu T., et al., 2016, ApJS, 222, 7Marleau G.-D., Klahr H., Kuiper R., Mordasini C., 2017, ApJ,836, 221Murillo N. M., Lai S.-P., Bruderer S., Harsono D., van DishoeckE. F., 2013, A&A, 560, A103Padoan P., Haugbølle T., Nordlund ˚A., 2012, ApJ, 759, L27Padoan P., Federrath C., Chabrier G., Evans II N. J., JohnstoneD., Jørgensen J. K., McKee C. F., Nordlund ˚A., 2014a, Pro-tostars and Planets VI, pp 77–100Padoan P., Haugbølle T., Nordlund ˚A., 2014b, ApJ, 797, 32Padoan P., Pan L., Haugbølle T., Nordlund ˚A., 2016, ApJ, 822,11Padoan P., Haugbølle T., Nordlund ˚A., Frimann S., 2017, ApJ,840, 48Paxton B., Bildsten L., Dotter A., Herwig F., Lesaffre P., TimmesF., 2011, ApJS, 192, 3Paxton B., et al., 2013, ApJS, 208, 4Paxton B., et al., 2015, ApJS, 220, 15Penston M. V., 1969, MNRAS, 144, 425Prialnik D., Livio M., 1985, MNRAS, 216, 37Prodanovi´c T., Steigman G., Fields B. D., 2010, MNRAS, 406,1108Rogers F. J., Nayfonov A., 2002, ApJ, 576, 1064Safron E. J., et al., 2015, ApJ, 800, L5Saumon D., Chabrier G., van Horn H. M., 1995, ApJS, 99, 713MNRAS , 1–20 (2017) Jensen & Haugbølle
Scholz A., Froebrich D., Wood K., 2013, MNRAS, 430, 2910Shu F. H., 1977, ApJ, 214, 488Stahler S. W., Shu F. H., Taam R. E., 1980a, ApJ, 241, 637Stahler S. W., Shu F. H., Taam R. E., 1980b, ApJ, 242, 226Stahler S. W., Shu F. H., Taam R. E., 1981, ApJ, 248, 727Teyssier R., 2002, A&A, 385, 337Tobin J. J., et al., 2016, ApJ, 818, 73Tosi M., 2000, in da Silva L., de Medeiros R., Spite M., eds, IAUSymposium Vol. 198, The Light Elements and their Evolution.p. 525Vasileiadis A., Nordlund ˚A., Bizzarro M., 2013, ApJ, 769, L8Vaytet N., Haugbølle T., 2017, A&A, 598, A116Vorobyov E., Elbakyan V., Hosokawa T., Sakurai Y., Guedel M.,Yorke H., 2017, preprint, ( arXiv:1706.00502 )Yen H.-W., Koch P. M., Takakuwa S., Krasnopolsky R., OhashiN., Aso Y., 2017, ApJ, 834, 178Zwintz K., 2015, Proceedings of the International AstronomicalUnion, 11, 552ˆa ˘A¸S559This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000