Evolution of circumstellar discs in young star-forming regions
Francisca Concha-Ramírez, Simon Portegies Zwart, Martijn J. C Wilhelm
MMNRAS , 1–12 (2020) Preprint 21 January 2021 Compiled using MNRAS L A TEX style file v3.0
Evolution of circumstellar discs in young star-forming regions
Francisca Concha-Ramírez ★ , Simon Portegies Zwart, Martijn J. C. Wilhelm Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
The evolution of circumstellar discs is highly influenced by their surroundings, in particular byexternal photoevaporation due to nearby stars and dynamical truncations. The impact of theseprocesses on disc populations depends on the dynamical evolution of the star-forming region.Here we implement a simple model of molecular cloud collapse and star formation to obtainprimordial positions and velocities of young stars and follow their evolution in time, includingthat of their circumstellar discs. Our disc model takes into account viscous evolution, internaland external photoevaporation, dust evolution, and dynamical truncations. The disc evolutionis resolved simultaneously with the star cluster dynamics and stellar evolution. Our resultsshow that an extended period of star formation allows for massive discs formed later in thesimulations to survive for several million years. This could explain massive discs surviving inregions of high UV radiation.
Key words: methods: numerical; planets and satellites: formation; stars: kinematics anddynamics;
Circumstellar discs are a natural consequence of the star formationprocess, and emerge within the first 10 yr after star formation(Williams & Cieza 2011). The star formation environment, rich ingas and newly-formed stars, can greatly affect the evolution of thediscs. The imprints that this environment will leave on the youngdiscs have important repercussions on their potential to form planetsand on the configurations of the planetary systems that eventuallydo develop.There are several ways in which the environment can influ-ence the evolution of circumstellar discs. Close encounters betweencircumstellar discs and stellar fly-bys can affect the size, mass, andsurface density of the discs. Close encounters can remove mass fromthe outskirts of the discs, decreasing both their mass and radius (e.g.Clarke & Pringle 1991, 1993; Pfalzner et al. 2005b; Breslau et al.2014; Vincke et al. 2015; Vincke & Pfalzner 2016; Portegies Zwart2016; Vincke & Pfalzner 2018; Cuello et al. 2018; Winter et al.2018a; Concha-Ramírez et al. 2019a). Several numerical imple-mentations of this process have shown that close encounters canlead to a hardening of the discs surface density (Rosotti et al. 2014),the formation of spiral arms and other structures (Pfalzner 2003;Pfalzner et al. 2005a), accretion bursts onto the host star (Pfalzneret al. 2008), and exchange of mass between discs (Pfalzner et al.2005b; Jílková et al. 2016). Observational evidence for the effectsof stellar fly-bys has been presented in several studies. Cabrit et al.(2006) study the ∼
600 au trailing “tail” in the disc of RW Aur Aand suggest that it might have been caused by a recent fly-by. Reche ★ E-mail: [email protected] et al. (2009) suggest that the spiral arms observed in the disc ofthe triple star system HD 141569 might be the result of a fly-by.Observations by Rodriguez et al. (2018) reveal newly-detected tidalstreams in RW Aur A and they propose that these might be the resultof many subsequent close encounters. Winter et al. (2018c) simu-late the disc around DO Tau, which presents a tidal tail, and arguethat this shape could have been caused by a close encounter withthe nearby triple system HV Tau. There is also evidence that theyoung disc of the solar system was affected by such an encounter.The sharp edge of the solar system at ∼
30 au could be a sign that apassing star truncated its early disc (Breslau et al. 2014; Punzo et al.2014). The highly eccentric and inclined orbits of the
Sednitos , agroup of 13 detected planetoids in the outskirts of the solar system,suggest they might have been captured from the disc of a passingnearby star (Jílková et al. 2015).Another mechanism that can alter the evolution of circum-stellar discs is photoevaporation. Photoevaporation is the processin which high energy photons heat the disc surface, causing themto evaporate. The source of these photons can be the host star(internal photoevaporation) or bright stars in the vicinity (exter-nal photoevaporation). Photoevaporation is driven by far ultraviolet(FUV), extreme ultraviolet (EUV), and X-ray photons (Johnstoneet al. 1998; Adams et al. 2004). The effects of internal and externalphotoevaporation on circumstellar discs are rather distinct. Internalphotoevaporation can clear areas of the disc at specific disc radii,causing the opening of gaps (Gorti et al. 2009; Gorti & Hollen-bach 2009; Owen et al. 2010; Font et al. 2004; Fatuzzo & Adams2008; Hollenbach et al. 2000). External photoevaporation can re-move mass from all over the disc surface, but the outer regions of © a r X i v : . [ a s t r o - ph . GA ] J a n Concha-Ramírez et al. the discs are more vulnerable because the material is less stronglybound to the host star (Johnstone et al. 1998; Adams et al. 2004).Observational evidence of external photoevaporation was firstobtained through the imaging of evaporating discs in the Orion neb-ula (O’dell & Wen 1994; O’dell 1998). These objects, now known as‘proplyds’, are circumstellar discs immersed in the radiation fieldsof nearby stars. Their cometary tail-like structure reveals the ongo-ing mass loss. Subsequent observations of the region showed thatcircumstellar disc masses decrease when close to massive stars. Thiseffect has been observed in several regions such as the Trapeziumcluster (e.g. Vicente & Alves 2005; Eisner & Carpenter 2006; Mannet al. 2014), the Orion Nebula Cluster (e.g. Mann & Williams 2010;Eisner et al. 2018), Cygnus OB2 (Guarcello et al. 2016), NGC 1977(Kim et al. 2016), NGC 2244 (Balog et al. 2007), Pismis 24 (Fanget al. 2012), NGC 2024 (van Terwisga et al. 2020), 𝜎 Orionis (Ans-dell et al. 2017), and 𝜆 Orionis (Ansdell et al. 2020). Younger andlow-mass star-forming regions such as Lupus, Taurus, Ophiuchus,and the Orion Molecular Cloud 2 tend to have higher average discmasses than denser regions such as the Orion Nebula Cluster (Eis-ner et al. 2008; Ansdell et al. 2016; Eisner et al. 2018; van Terwisgaet al. 2019). van Terwisga et al. (2020) present the discovery oftwo distinct disc populations, in terms of mass, in the NGC 2024region. The discs to the east of the region are embedded in a densemolecular ridge and are more massive than the discs outside theridge, which are also closer to two OB type stars. They propose thatthe difference in masses is caused by the eastern population beingprotected from the radiation of nearby massive star IRS 1.Several models have demonstrated external photoevaporationis efficient in depleting disc masses on timescales much shorterthan their estimated lifetimes of ∼
10 Myr (e.g. Scally & Clarke2001; Adams et al. 2006; Fatuzzo & Adams 2008; Haworth et al.2016), even in low radiation fields (Facchini et al. 2016; Kim et al.2016; Haworth et al. 2017). Because external photoevaporation iscaused by massive stars in the vicinity of the discs, the extent of itseffects depends on the density of the stellar region and the numberof massive stars in the surroundings. Even in high density regions(N ∗ (cid:38) pc − ) the disc mass-loss rates caused by external pho-toevaporation are orders of magnitude higher than those causedby dynamical truncations (Winter et al. 2018b, 2019b; Concha-Ramírez et al. 2019b). Concha-Ramírez et al. (2021) show that, inregions of local stellar densities N ∗ >
100 pc − , external photoe-vaporation can evaporate up to 90% of circumstellar discs within2 . ∼
10 M (cid:12) pc − ), only ∼
60% ofdiscs are evaporated within the same timescale. Winter et al. (2020)model a region comparable to the central molecular zone of theMilky Way (surface density Σ = M (cid:12) pc − ) and find that ex-ternal photoevaporation destroys 90% of circumstellar discs within1 . Σ =
12 M (cid:12) pc − ) they find amean disc dispersion timescale of 3 . ? who find that external photoevaporation destroys 50%of discs within 1 . ∼
100 M (cid:12) pc − , andwithin 2 . ∼
10 M (cid:12) pc − .While observational and numerical evidence indicate that discmasses decrease with increasing stellar density, massive discs arestill observed in high density regions. In particular, there are discsin the ONC whose masses are much higher than other discs in theproximity of the massive star 𝜃 Ori C. If the discs are coeval with 𝜃 Ori C they should have already been dispersed by external phote-vaporation, unless they were extraordinarily massive to begin with(M disc (cid:38) (cid:12) ). Alternatively, 𝜃 Ori C would have to be consid-erably younger than the ONC average ( (cid:46) . We model several different astrophysical processes which operatesimultaneously and at very distinct scales: the collapse of a molecu-lar cloud, star formation, stellar dynamics, and viscous circumstellardiscs which are affected by dynamical truncations and photevapora-tion. We bring these processes together using the Astrophysical Mul-tipurpose Software Environment, AMUSE (Portegies Zwart et al.2013; Pelupessy et al. 2013). The results presented in this workare obtained through two different simulation stages: first, we per-form a simple model of the collapse of a molecular cloud and thesubsequent star formation process. This returns a spatial, velocity,and mass distribution of stars to be used in the second simulationstage, in which we follow the evolution of the circumstellar discs in
MNRAS000
MNRAS000 , 1–12 (2020) ircumstellar discs in young star-forming regions the star-forming regions. This second stage encompasses the stel-lar dynamics, stellar evolution, viscous evolution of the discs, andphotoevaporation. All the code developed for this work is availableonline . The first stage of the simulations deals with a simple model of thestar formation process. We simulate the collapse of a molecularcloud using the smoothed particle hydrodynamics (SPH) code Fi (Pelupessy et al. 2004). We model the star formation process throughthe use of sink particles, which are created from regions of the cloudwhere the local gas density is higher than a threshold. We set thisthreshold at 1M (cid:12) / 𝜖 , where 𝜖 = .
05 pc is the softening of thesimulation. Once a sink particle forms it continues accreting gasfrom the molecular cloud. Each of these sink particles can formseveral stars.Since we look to preserve a power-law stellar mass functionsimilar to the one observed in the galaxy, the stars in our simulationscannot be formed without introducing a predetermined star forma-tion efficiency (SFE). We set a SFE of 0.3 (Lada & Lada 2003). Weimplement this SFE by keeping track of the mass in sinks, since thisis the mass that will eventually be turned into stars. When the totalsinks mass reaches 30% of the mass of the cloud, the SPH code isstopped. We keep track of the existing sinks in a separate dynamicscode to guarantee that they continue to move after the SPH code isstopped. The star formation process continues until all the mass inthe sinks has turned into stars.The star formation process begins once sink particles haveaccreted enough mass to sample stars from a random initial massfunction (IMF). We base the star formation mechanism on Wallet al. (2019). We begin by sampling a random stellar mass 𝑚 froma Kroupa IMF (Kroupa 2001) of 10.000 stars, with lower limit 0.08M (cid:12) and upper limit 150 M (cid:12) (Wall et al. 2019). Then we check ifthere is a sink massive enough to form a star of mass 𝑚 . If there isone, we subtract the mass 𝑚 from the sink and create a star particle.If 𝑚 ≤ . (cid:12) the star will have a circumstellar disc (see section2.2.2), and we subtract the mass of the star and the initial mass ofthe disc from the sink. The position of the newly formed star isdetermined by taking the position of the sink and adding a randomoffset in each spatial dimension. This offset is calculated within thesink diameter. The velocity of the new star is set to the velocity ofits birth sink.After a sink has formed a star, we set a delay time whichmust pass before it creates a new star. We implement this step tocounteract the fact that our sampling will be biased toward forminglow mass stars, by allowing sinks to become more massive beforeforming another star. This delay is implemented as an exponentiallydecaying timescale:t delay = t ff exp (cid:18) − . (cid:19) (1)where t ff is the free-fall time scale of the corresponding sink and 𝑡 is the current model time.During the molecular cloud collapse simulation we keep trackof the mass, position, velocity, and birth time of the newly formedstars. This data will then be given as input for the second part of thesimulations, which involve the disc evolution and stellar dynamics. https://doi.org/10.5281/zenodo.4436316 The dynamical evolution of the stars is calculated using the 4 th -orderHermite integrator ph4 , which is evolved simultaneously with theSPH integrator in a leap-frog scheme using the Bridge (Fujii et al.2007) coupling method in AMUSE (see Portegies Zwart et al. 2020,for implementation details).
The second simulation stage begins at the time when the first starhas formed. For each star, we evolve its disc and calculate its exter-nal photoevaporation mass loss rate as explained in sections 2.2.1and 2.2.2, respectively. The star formation process ends when all themass in sinks has been formed into stars. Then, the leftover gas isexpelled and we only deal with the stellar dynamics and processesexplained in the following sections. This second stage of the sim-ulations evolve for 2 Myr after the last star has formed. Below wedescribe how each of the processes is modeled.
We model circumstellar discs using the Viscous Accretion discEvolution Resource (VADER, Krumholz & Forbes 2015). VADERmodels the viscous transport of mass and angular momentum on athin, axisymmetric disc. We use the standard disc profile of Lynden-Bell & Pringle (1974) to establish the initial column density of thediscs as: Σ ( 𝑟, 𝑡 = ) ≈ 𝑚 𝑑 𝜋𝑟 𝑑 (cid:0) − 𝑒 − (cid:1) exp (− 𝑟 / 𝑟 𝑑 ) 𝑟 (2)where 𝑟 𝑑 is the initial disc radius, 𝑚 𝑑 is the initial disc mass, and Σ is a normalisation constant. We consider the characteristic radius tobe 𝑟 𝑐 ≈ 𝑟 𝑑 (Anderson et al. 2013).For the external photoevaporation process we keep track ofthe outer edge of the disc. We define the disc radius as the radiuswhich encloses 99% of the disc mass (Anderson et al. 2013). Themass loss due to external photoevaporation (section 2.2.2), as wellas dynamical truncations (section 2.2.5), causes the disc to developa steep density profile at the outer edge. The location of the edgeis insensitive to the value of Σ edge , given that it is sufficiently low(Clarke 2007). We set the column density outside 𝑟 𝑑 to a negligiblevalue Σ edge = − g cm − .Each of the VADER discs in the simulations consists of a gridof 100 logarithmically spaced cells, ranging from 0 .
05 au to 2000 au.All the discs have a constant turbulence parameter of 𝛼 = × − . We calculate the mass loss due to external photoevaporation usingthe Far-ultraviolet Radiation Induced Evaporation of Discs (FRIED)grid (Haworth et al. 2018b). This grid provides a pre-calculated setof mass loss rates for discs immersed in UV radiation fields ofvarying strengths, from 10 G to 10 G , where G is the FUVfield measured in Habing units, 1 . × − erg s − cm − (Habing1968). The grid spans discs of mass ∼ − M Jup to 10 M Jup ,radius from 1 au to 400 au, and host star mass from 0 .
05 M (cid:12) to1 . (cid:12) .To stay inside the boundaries of the FRIED grid we consideronly stars of mass M ∗ ≤ . (cid:12) as having a circumstellar disc.More massive stars are considered to generate radiation. This massdistinction is for external photoevaporation calculations only; for MNRAS , 1–12 (2020)
Concha-Ramírez et al. the stellar dynamics evolution there is no separation between thesetwo stellar groups.Far-ultraviolet (FUV) photons dominate in the external photo-evaporation process (Armitage 2000; Adams et al. 2004; Gorti &Hollenbach 2009). We calculate the FUV radiation from the mas-sive stars by pre-computing a relation between stellar mass andFUV luminosity using the UVBLUE spectral library (Rodriguez-Merino et al. 2005). The obtained fit is presented in Figure 2 ofConcha-Ramírez et al. (2019b). We use that fit to determine theFUV radiation emitted by each massive star at every simulationtime step.We model the external photoevaporation process as follows.At every time step we calculate the distance from every disc toevery star of mass M ∗ > . (cid:12) and determine the total radia-tion received by each disc. We do not consider extinction due tointerstellar material. We interpolate from the FRIED grid using thecalculated total radiation and the disc parameters to find a photoe-vaporation driven mass loss rate for each disc. Assuming the massloss rate (cid:164) M to be constant during the current time step, we use it tocalculate the total mass loss. This mass is then removed from theouter regions of the disc. We move through the disc cells startingfrom the outermost one, removing mass from each, until the requiredamount of mass has been removed. External photoevaporation thenresults in a decrease of disc mass and disc radius.Under certain circumstances, external photoevaporation canbe dominated by extreme ultraviolet (EUV) photons. This happenswhen a disc is closer to a radiating star than a minimum distancederived by Johnstone et al. (1998) as: 𝑑 𝑚𝑖𝑛 (cid:39) × (cid:18) 𝜀 𝑓 𝑟 Φ (cid:19) − / r / 𝑑 cm (3)where 𝑓 𝑟 is the fraction of EUV photons absorbed in the ionizingflow, Φ = Φ 𝑖 s − is the EUV luminosity of the source, 𝜀 is adimensionless normalizing parameter, (cid:16) 𝜀 𝑓 𝑟 Φ (cid:17) / ≈
4, and 𝑟 𝑑 = 𝑟 𝑑 𝑐𝑚 with 𝑟 𝑑 the disc radius.When the distance 𝑑 between a disc and a massive star is 𝑑 < 𝑑 𝑚𝑖𝑛 , the disc enters the EUV-dominated photoevaporationregime. The mass loss in this case is calculated as: (cid:164) 𝑀 𝐸𝑈𝑉 = . × − ( + 𝑥 ) 𝑥 𝜀𝑟 𝑑 M (cid:12) yr − (4)with 𝑥 ≈ . 𝜀 ≈ / cm (Pascucci et al. 2016). After a disc is dispersed, its hoststar continues to evolve normally in the stellar dynamics code. Internal photoevaporation is driven by X-ray radiation (Owen et al.2010, 2012). We calculate the X-ray luminosity of the stars withdiscs using the fit obtained by Flaccomio et al. (2012) for T-Tauristars as a function of stellar mass:log (cid:18) 𝐿 𝑋 erg s − (cid:19) = . (cid:18) M ∗ (cid:12) (cid:19) +
30 (5)where M ∗ is the mass of the host star.Picogna et al. (2019) calculate X-ray mass loss profiles and mass loss rates for a star of mass 0.7 M (cid:12) . Owen et al. (2012)developed scaling relations that allow to calculate these values forstars with masses M ∗ ≤ . (cid:12) . We combine the results of Picognaet al. (2019) with the scaling relations of Owen et al. (2012) to spana larger range of stellar masses. The internal photoevaporation massloss rate is then given by: (cid:164) 𝑀 𝑋 = Δ (cid:18) M ∗ . (cid:12) (cid:19) − . M (cid:12) yr − , (6)where Δ = − . (cid:20) ( ln ( log ( 𝐿 𝑋 )) − . ) . × − (cid:21) − . (cid:164) Σ 𝑤 ( 𝑥 ) = ln ( ) (cid:32) 𝑎 ln ( 𝑥 ) 𝑥 ln ( ) 𝑏 + 𝑏 ln ( 𝑥 ) 𝑥 ln ( ) + 𝑐 ln ( 𝑥 ) 𝑥 ln ( ) + 𝑑 ln ( 𝑥 ) 𝑥 ln ( ) + 𝑒 ln ( 𝑥 ) 𝑥 ln ( ) + 𝑓𝑥 ln ( ) (cid:33) (cid:164) 𝑀 𝑤 ( 𝑥 ) 𝜋𝑥 M (cid:12) au − yr − , (8)where (cid:164) 𝑀 𝑤 ( 𝑥 ) = (cid:164) 𝑀 𝑋 𝑎 log 𝑥 + 𝑏 log 𝑥 + 𝑐 log 𝑥 + 𝑑 log 𝑥 + 𝑒 log 𝑥 + 𝑓 log 𝑥 + 𝑔 (9)with 𝑎 = − . 𝑏 = . 𝑐 = − . 𝑑 = . 𝑒 = − . 𝑓 = . 𝑔 = − . 𝑥 = . (cid:18) 𝑟 au (cid:19) (cid:18) M ∗ (cid:12) (cid:19) − (10)is the scaling from Owen et al. (2012), where M ∗ is the mass of thehost star.The internal photoevaporation process removes mass from thedisc following the profile defined in Eq. 8. In the case where agrid cell contains less mass than is prescribed to be removed, thisexcess is removed in the nearest outer cell. As the cells are traversedinside-out, this takes the form of inside-out disc clearing. Circumstellar discs are composed of gas and dust. Initially, a 100:1gas:dust ratio is assumed for the composition of the discs. This valueis derived from the consideration that the ratio is inherited from theinterstellar medium (Bohlin et al. 1978). Grain growth might resultin much lower gas:dust ratios (Williams & Best 2014), and the ratiois likely to change during the lifetime of a disc (Manara et al. 2020).Models of dust evolution and radial drift (Birnstiel et al. 2010;Rosotti et al. 2019) show that the dust:gas ratio decreases in time.In the present work we introduce a simple prescription for the dustevolution inside circumstellar discs.We follow the prescription of Haworth et al. (2018a) to cal-culate the mass loss rate of dust entrapped in the photoevaporationwind. This mass loss rate is described as: (cid:164) 𝑀 𝑑𝑤 = 𝛿 (cid:164) 𝑀 / 𝑔𝑎𝑠 (cid:18) 𝜈 𝑡ℎ 𝜋𝐹𝐺 M ∗ 𝜌 𝑔 𝑎 𝑚𝑖𝑛 (cid:19) / exp (cid:32) − 𝛿 ( 𝐺 M ) / 𝑡 𝑑 / (cid:33) , (11) MNRAS000
30 (5)where M ∗ is the mass of the host star.Picogna et al. (2019) calculate X-ray mass loss profiles and mass loss rates for a star of mass 0.7 M (cid:12) . Owen et al. (2012)developed scaling relations that allow to calculate these values forstars with masses M ∗ ≤ . (cid:12) . We combine the results of Picognaet al. (2019) with the scaling relations of Owen et al. (2012) to spana larger range of stellar masses. The internal photoevaporation massloss rate is then given by: (cid:164) 𝑀 𝑋 = Δ (cid:18) M ∗ . (cid:12) (cid:19) − . M (cid:12) yr − , (6)where Δ = − . (cid:20) ( ln ( log ( 𝐿 𝑋 )) − . ) . × − (cid:21) − . (cid:164) Σ 𝑤 ( 𝑥 ) = ln ( ) (cid:32) 𝑎 ln ( 𝑥 ) 𝑥 ln ( ) 𝑏 + 𝑏 ln ( 𝑥 ) 𝑥 ln ( ) + 𝑐 ln ( 𝑥 ) 𝑥 ln ( ) + 𝑑 ln ( 𝑥 ) 𝑥 ln ( ) + 𝑒 ln ( 𝑥 ) 𝑥 ln ( ) + 𝑓𝑥 ln ( ) (cid:33) (cid:164) 𝑀 𝑤 ( 𝑥 ) 𝜋𝑥 M (cid:12) au − yr − , (8)where (cid:164) 𝑀 𝑤 ( 𝑥 ) = (cid:164) 𝑀 𝑋 𝑎 log 𝑥 + 𝑏 log 𝑥 + 𝑐 log 𝑥 + 𝑑 log 𝑥 + 𝑒 log 𝑥 + 𝑓 log 𝑥 + 𝑔 (9)with 𝑎 = − . 𝑏 = . 𝑐 = − . 𝑑 = . 𝑒 = − . 𝑓 = . 𝑔 = − . 𝑥 = . (cid:18) 𝑟 au (cid:19) (cid:18) M ∗ (cid:12) (cid:19) − (10)is the scaling from Owen et al. (2012), where M ∗ is the mass of thehost star.The internal photoevaporation process removes mass from thedisc following the profile defined in Eq. 8. In the case where agrid cell contains less mass than is prescribed to be removed, thisexcess is removed in the nearest outer cell. As the cells are traversedinside-out, this takes the form of inside-out disc clearing. Circumstellar discs are composed of gas and dust. Initially, a 100:1gas:dust ratio is assumed for the composition of the discs. This valueis derived from the consideration that the ratio is inherited from theinterstellar medium (Bohlin et al. 1978). Grain growth might resultin much lower gas:dust ratios (Williams & Best 2014), and the ratiois likely to change during the lifetime of a disc (Manara et al. 2020).Models of dust evolution and radial drift (Birnstiel et al. 2010;Rosotti et al. 2019) show that the dust:gas ratio decreases in time.In the present work we introduce a simple prescription for the dustevolution inside circumstellar discs.We follow the prescription of Haworth et al. (2018a) to cal-culate the mass loss rate of dust entrapped in the photoevaporationwind. This mass loss rate is described as: (cid:164) 𝑀 𝑑𝑤 = 𝛿 (cid:164) 𝑀 / 𝑔𝑎𝑠 (cid:18) 𝜈 𝑡ℎ 𝜋𝐹𝐺 M ∗ 𝜌 𝑔 𝑎 𝑚𝑖𝑛 (cid:19) / exp (cid:32) − 𝛿 ( 𝐺 M ) / 𝑡 𝑑 / (cid:33) , (11) MNRAS000 , 1–12 (2020) ircumstellar discs in young star-forming regions where 𝛿 is the initial dust:gas ratio (10 − ), (cid:164) 𝑀 𝑔𝑎𝑠 is the gas massloss rate (determined as explained in sections 2.2.2 and 2.2.3), 𝜈 𝑡ℎ = √︁ 𝑘 𝑏 𝑇 / 𝜋 / 𝜇 / 𝑚 𝐻 is the mean thermal speed of the gas particles, 𝐹 is the solid angle subtended by the disc at the outer edge, 𝜌 𝑔 is the grain mass density (1 g / cm , Facchini et al. (2016)), and 𝑎 𝑚𝑖𝑛 is the minimum grain size at the disc radius R 𝑑 . We assume 𝑎 𝑚𝑖𝑛 = . 𝜇 m (Haworth et al. 2018a; Facchini et al. 2016).This model takes into account what fraction of the dust isentrained in the photoevaporation wind, and how it decreases overtime due to dust growth. Mass is removed from a single, scalarreservoir. The radial structure is implicitly assumed to follow the gasstructure, multiplied by the dust-to-gas ratio 𝛿 ∼ .
01, and doesn’taccount for the dust fraction enhancement due to the evaporation-resistant dust population.
We calculate a semi-analytical truncation radius based on Adams(2010), who propose that the new radius of a disc after a truncatingencounter is 𝑅 (cid:48) ≈ 𝑏 / 𝑏 is the pericentre distance of theencounter. We combine this with the mass dependence of Breslauet al. (2014) to define a truncation radius: 𝑅 (cid:48) = 𝑟 𝑒𝑛𝑐 (cid:18) 𝑚 𝑚 (cid:19) . , (12)where 𝑚 and 𝑚 are the masses of the encountering stars. We follow(Portegies Zwart 2016) in defining an intial collisional radius of 𝑟 𝑐𝑜𝑙 = .
02 pc for all stars. This value is updated to 𝑟 𝑐𝑜𝑙 = . 𝑟 𝑒𝑛𝑐 after every encounter, to make sure that each encounter is onlydetected once within the time step. Not all encounters result in disctruncation. If the calculated truncation radius 𝑅 (cid:48) is larger than thecurrent radius of an encountering disc, the disc is not affected by theencounter. If a disc is truncated in an encounter, we set the new radiusof a disc to 𝑅 (cid:48) by making the column density Σ edge = − g cm − for every disc cell outside 𝑅 (cid:48) . The truncated disc then continues toexpand viscously.Dynamical encounters not only change the disc sizes and stripmass from the outskirts, but can also lead to changes in the massdistribution of the discs and mass exchange can occur between theencountering discs (Pfalzner et al. 2005b; Rosotti et al. 2014; Jílkováet al. 2015; Portegies Zwart 2016). Because of our implementationof the discs (section 2.2.1) we do not consider mass exchange or anyother changes to the mass distribution during dynamical encounters,other than truncation. When a disc is truncated in our model, all themass outside its new radius 𝑅 (cid:48) is simply lost. Our simulations start with an spherical cloud model of mass 10 M (cid:12) and initial radius 3 pc. We use 32.000 SPH particles, whichresults in a resolution of 0.3 M (cid:12) per particle. The softening in thesimulations is 0 .
05 pc. We use a power-law velocity spectrum tomodel large scale turbulence (Bate et al. 2003). Each realization ofa molecular cloud has a different random seed, so the substructurethat originates is different in every run of the simulation. We run6 realizations of the molecular cloud collapse simulations. Thererealizations differ only in the random seed used to determine theposition and velocities of the SPH particles.
We choose the initial disc radii as: 𝑟 𝑑 ( 𝑡 = ) = 𝑅 (cid:48) (cid:18) 𝑀 ∗ 𝑀 (cid:12) (cid:19) . , (13)where 𝑅 (cid:48) is a constant. We choose 𝑅 (cid:48) =
30 au, which results ininitial disc radii between ∼ ∼
40 au. This is in agreementwith observations that suggest that young circumstellar discs aregenerally quite compact (radii around 20 to 50 au, (Trapman et al.2020; Tobin et al. 2020)). We set the initial radius of each disc onthe grid through the procedure explained in section 2.2.1.The initial mass of the discs is:M 𝑑 ( 𝑡 = ) = . ∗ . (14)This yields initial disc masses ranging from ∼ 𝐽𝑢 𝑝 to ∼
200 M
𝐽𝑢 𝑝 . In Figure 2 we show the number of stars in time for each simulationrun. In Table 1 we present the final number of stars in each run.The mean number of stars created in six runs is 6146 ± 𝑄 = 𝑚𝑠 , (15)where 𝑚 is the mean length of the minimum spanning tree and 𝑠 is the mean separation between the stars. Regions with Q > . < . 𝑀 ∗ > . (cid:12) . In Table 2 we summarize the values for Q and theestimated ages for each region, along with the corresponding refer-ences. The simulation results span a range of Q ∼ . ∼ . MNRAS , 1–12 (2020)
Concha-Ramírez et al.
Figure 1.
Evolution of the molecular cloud collapse and star formation process in run
Figure 2.
Number of stars in time for each simulation run.
Figure 3.
Virial radius of the simulations in time. The solid lines correspondto the virial radius while star formation is still ongoing, and the dotted linesafter it has ended. The black line shows the virial radius in time of a Plummersphere with 6000 stars and initial virial radius 3 pc for comparison.
Figure 4.
Q parameter of our simulations in time, and values for observedstar-forming regions. The Q parameter in our simulations considers onlystars with masses 𝑀 ∗ > . (cid:12) . The solid lines represent the times wherestar formation is still ongoing, and the dotted lines after all stars have beenformed. and Taurus (Cartwright & Whitworth 2004) are highly substruc-tured. In observations of star-forming regions, Q might vary de-pending on membership uncertainty (Parker & Meyer 2012). Thisleads to some regions having more than one value of Q, such asCorona Australis, Chamaeleon, the ONC, Ophiuchus, and UpperScorpio.As can be seen in Figure 4, the final shapes of the clustersgenerated by our simulations are smoother than those of observedclusters. This might be caused by the fact that we do not take intoaccount stellar winds and feedback from stellar processes, whichmight be important for the emergence of substructure (e.g. Mac Low& Klessen 2004; Hansen et al. 2012; Offner & Arce 2015). Weexpand on this discussion in Section 4.Another way to quantify the structure of a star-forming regionis by measuring its fractal dimension, 𝐹 𝑑 . In Figure 5 we show theevolution of the fractal dimension in time for each simulation run,along with a Plummer sphere of 6000 stars and initial virial radius3 pc for comparison. The solid lines show the fractal dimensionwhile star formation is still ongoing. The dotted lines show the MNRAS000
Q parameter of our simulations in time, and values for observedstar-forming regions. The Q parameter in our simulations considers onlystars with masses 𝑀 ∗ > . (cid:12) . The solid lines represent the times wherestar formation is still ongoing, and the dotted lines after all stars have beenformed. and Taurus (Cartwright & Whitworth 2004) are highly substruc-tured. In observations of star-forming regions, Q might vary de-pending on membership uncertainty (Parker & Meyer 2012). Thisleads to some regions having more than one value of Q, such asCorona Australis, Chamaeleon, the ONC, Ophiuchus, and UpperScorpio.As can be seen in Figure 4, the final shapes of the clustersgenerated by our simulations are smoother than those of observedclusters. This might be caused by the fact that we do not take intoaccount stellar winds and feedback from stellar processes, whichmight be important for the emergence of substructure (e.g. Mac Low& Klessen 2004; Hansen et al. 2012; Offner & Arce 2015). Weexpand on this discussion in Section 4.Another way to quantify the structure of a star-forming regionis by measuring its fractal dimension, 𝐹 𝑑 . In Figure 5 we show theevolution of the fractal dimension in time for each simulation run,along with a Plummer sphere of 6000 stars and initial virial radius3 pc for comparison. The solid lines show the fractal dimensionwhile star formation is still ongoing. The dotted lines show the MNRAS000 , 1–12 (2020) ircumstellar discs in young star-forming regions Table 1.
Run number, final number of stars ( 𝑁 ∗ ), time of the end of star for-mation (t SFend ), Q parameter, and fractal dimension ( 𝐹 _ 𝑑 ) for our simulationresults. Region 𝑁 ∗ t SFend [Myr] Q end 𝐹 𝑑 Run
Figure 5.
Fractal dimension of the simulations in time. The solid linescorrespond to the fractal dimension while star formation is still ongoing,and the dotted lines after it has ended. The black line shows the fractaldimension in time of a Plummer sphere with 6000 stars and initial virialradius 3 pc for comparison. evolution after the star formation process has ended. Once a regionreaches a stable value of 𝐹 𝑑 , it remains the same up to the endof the simulation. We also show measured values of the fractaldimension for several observed regions. As with the Q parameter,when measuring the fractal dimension from observations the valuesmight vary depending on membership, leading to regions with morethan one value of 𝐹 𝑑 . In Table 2 we present the values of 𝐹 𝑑 for allobserved regions.From Figures 3, 4, and 5 it can be seen that after the gas in theregions is expelled, the In Figure 6 we show the disc mass versus local number density forall simulation runs. The left panel shows the discs at the end ofthe star formation process. The right panel shows the discs at theend of the simulations, 2 Myr after star formation has finished. Thecolor of each disc represents the time at which it was born. Wecalculate the local stellar density using the method by Casertano &Hut (1985) with the five nearest neighbours. While star formationis still ongoing, stars and discs form in regions spanning the wholerange of stellar density. In particular, given our implementation ofthe star formation process with sink particles, many stars tend toform in regions of high stellar densities. Once the star formationprocess ends and the gas is expelled, the clusters expand to regainvirial equilibrium. This brings about an overall decrease in localnumber density. As stars are no longer being formed, the areas of high stellar density become less populated, since discs are losingmass due to photoevaporation. By the end of the simulations, almostno discs are present in regions of local density (cid:38) stars pc − .Even in low density regions, the maximum disc mass decreases from ∼
200 M
𝐽𝑢 𝑝 to ∼
100 M
𝐽𝑢 𝑝 . Massive, large discs are present onlyin regions of low stellar density. The vertical ridges seen in the leftpanel of Figure 6 are a numerical effect on the measurement of discmass and do not affect the overall results.In Figure 7 we show the mean mass in time for the discs ineach run. The solid lines show the times while star formation is stillongoing, and the dotted lines once the last star has formed. Whilestar formation is still happening, there is variability in the meandisc mass given the fact that some discs are exposed to externalphotoevaporation and some new discs are being formed. However,once star formation stops, the decrease in disc mass is sharp.In Figure 8 we present the binned mean gas and dust masses ofthe discs versus the projected local stellar number density, at the endof the simulations. The local stellar density is calculated in the sameway as for Figure 6, but with the distances between stars projectedto two dimensions. The binned mean is calculated using a rollingbin spanning 100 stars. The dust mass remains relatively constantacross densities, whereas the gas mass shows a slight decrease withstellar density. This can be explained by the fact that some dust islost in photoevaporation early on, but it soon becomes resilient tothis effect. Gas, on the other hand, is consistently lost throughoutthe simulations.In Figure 9 we show the cumulative distribution of disc dustmasses at the end of star formation and at the end of the simulations.The lines show the mean values for all runs and the shaded areasshow the standard deviation. By the end of star formation, all existingdiscs have dust masses (cid:38) ⊕ , up to ∼ ⊕ . After 2 Myrof evolution, the distribution spans discs with dust masses from ∼ − M ⊕ to ∼
400 M ⊕ . Around 80% of the final discs havedust masses higher than 10 M ⊕ , which is a lower limit mass for theformation of rocky planets and the cores of gas giants (Ansdell et al.2016). The masses obtained in the present simulations are higherthan those in Concha-Ramírez et al. (2021). This suggests that theextended period of star formation might play a role in allowingmassive discs to survive for longer. We expand on this discussionin Section 4.In Figure 10 we show the evolution of the mean dust-to-gasdisc mass ratio, or 𝛿 = M dust / M gas in time, aggregated for all simu-lation runs. By definition, the discs in the simulations are initializedwith a value of 𝛿 = − . As gas mass is lost from the discs dueto photoevaporation, this ratio tends to increase in time. By theend of the simulations the mean value is 𝛿 ≈ × − . The dustquickly becomes resistant to photoevaporation, but gas continuesbeing removed from the discs. In previous work (Concha-Ramírez et al. 2019b, 2021) we per-formed simulations of star clusters where stars with masses M ∗ ≤ . (cid:12) have circumstellar discs. These discs were subject to vis-cous evolution and external photoevaporation. We assumed the dustmass to be 1% of the total disc mass, and that only gas mass isremoved through photoevaporation. In those previous simulationswe found that discs get depleted of mass quickly enough so that ∼
60% to 90% of them (depending on the stellar density of theregion) are completely dispersed within 2 Myr of evolution. Theinitial spatial distribution of the stars were Plummer spheres, and
MNRAS , 1–12 (2020)
Concha-Ramírez et al.
Table 2.
Region name, number of stars ( 𝑁 ∗ ), age, Q parameter, and fractal dimension ( 𝐹 _ 𝑑 ) for our simulation results and observed regions. References: (a)Parker (2014); (b) Neuhäuser & Forbrich (2008); (c) Luhman et al. (2016); (d) Parker & Alves de Oliveira (2017); (e) Hillenbrand & Hartmann (1998); (f)Cartwright & Whitworth (2004); (g) Hartmann (2002); (h) Kraus & Hillenbrand (2008); (i) Simon (1997); (j) Bontemps et al. (2001); (k) Luhman (2007); (l)Wright et al. (2010); (m) Wright et al. (2014); (n) Carpenter et al. (2006); (o) Luhman & Mamajek (2012); (p) Sacco et al. (2017); (q) Galli et al. (2020); (r)Luhman & Esplin (2020); (s) Luhman (2018). Region 𝑁 ∗ Age [Myr] Q 𝐹 𝑑 Corona Australis (CrA) ∼ ( 𝑞 ) ∼ . ( 𝑎 ) . ( 𝑏 ) , 0 . ( 𝑎 ) -NGC 1333 ∼ ( 𝑐 ) ∼ . ( 𝑐 ) . ( 𝑑 ) -ONC ∼ ( 𝑒 ) ∼ . ( 𝑒 ) . ( 𝑒 ) , 0 . ( 𝑎 ) -Taurus ∼ ( 𝑠 ) ∼ . ( 𝑓 ) . ( 𝑓 ) . ( 𝑓 ) , 1 . ± . ( 𝑔 ) , 1 . ± . ( ℎ ) , 1 . ± . ( 𝑖 ) Trapezium ∼ ( 𝑒 ) ∼ . ( 𝑒 ) - 1 . ± . ( 𝑖 ) Ophiuchus 199 ( 𝑓 ) . ± . ( 𝑗 ) . ( 𝑓 ) , 0 . ( 𝑎 ) . ± . ( 𝑖 ) Chamaeleon I 120 ( 𝑝 ) . ± . ( 𝑘 ) . ( 𝑓 ) , 0 . ( 𝑎 ) , 0 . ± . ( 𝑝 ) . ( 𝑓 ) Cygnus OB2 ∼ ( 𝑙 ) . ± . ( 𝑙 ) . ± . ( 𝑚 ) -IC 348 ∼ ( 𝑐 ) . ± . ( 𝑐 ) . ( 𝑑 ) -Upper Scorpio ∼ ( 𝑟 ) . ± . ( 𝑛 ) . ( ℎ ) , 0 . ( 𝑎 ) . ± . ( ℎ ) Figure 6.
Disc mass versus local number density for all discs in all simulation runs. The left panel shows the discs at the end of the star formation process. Theright panel shows the discs at the end of the simulations, 2 Myr after star formation has finished. The size of the points is proportional to the disc radius. Thecolor shows the time at which each disc was born. all stars were formed at the same time and evolved for the 2 Myrthat the simulations lasted.In the present work we improve our previous model in twomain ways. First, we implement a simple model of molecular cloudcollapse and star formation. This results in the stars not havingspherical initial distributions, and in new stars being added to thesimulation in time throughout an extended star formation process.Second, we implement a simple model for the evolution of the dustin the discs (Haworth et al. 2018a), in order to follow its progressiondirectly and not simply approximate it to 1% of the total mass. Wealso implement internal photoevaporation as an additional means toremove gas mass from the discs.During the first few million years of evolution the star formationprocess is ongoing and, as discs lose mass due to photoevaporation,the new discs that are constantly being formed keep the mean massof the overall population relatively unaffected (Figure 7). As soonas this constant replenishing of stars and discs stops, the mean mass of the discs quickly drops. This is consistent with the results fromConcha-Ramírez et al. (2019b) and Concha-Ramírez et al. (2021).However, disc masses at the end of the simulations (see Figure 8and the bottom panel of Figure 6) are higher than in our previoussimulations. In particular, unlike in our previous calculations, thedust masses in this work do not diminish as strongly with increasingstellar density.This discrepancy can be explained by an interplay of differentprocesses. The dust model that we implement is such that somedust is initially trapped in the photoevaporation wind. As the discscontinue to evolve, the dust becomes resistant to photoevaporationand it is mostly gas mass that is lost from the discs. After a fewhundred thousand years of evolution, the amount of dust mass be-comes independent of the radiation received by the discs. The gasmass does decrease with increasing stellar density (see top panel ofFigure 8) and in time (Figure 7). Ansdell et al. (2016) propose thata rapid depletion of gas mass in discs can lead to the observed char-
MNRAS000
MNRAS000 , 1–12 (2020) ircumstellar discs in young star-forming regions Figure 7.
Mean disc mass in time for each simulation run. The solid linescorrespond to ongoing star formation, and the dotted lines after it has ended.The standard deviation in two of the runs is shown as an example; the otherruns have standard deviations of similar magnitude.
Figure 8.
Binned mean gas mass (top panel, M
𝐽𝑢𝑝 ) and dust mass (bottompanel, M ⊕ ) at the end of the simulations versus projected local stellar density.The binned mean is calculated using a rolling bin spanning 100 stars. Thelocal stellar number density is calculated as specified in section 3.2. Figure 9.
Cumulative distribution of disc dust and gas masses at the end ofstar formation (dashed lines) and at the end of the simulations (solid lines).The lines show the mean values for all runs and the shaded areas show thestandard deviation.
Figure 10.
Mean M dust / M gas versus time. acteristics of exoplanet populations. Traditional theories of planetformation suggest that ∼ ⊕ cores should be able to rapidlyaccrete gaseous envelopes, and reach masses of ∼ 𝐽𝑢 𝑝 within ∼ . ∼ ⊕ cores form in discs which are already depletedof gas, this would result in the observed discrepancy in the numberof planetary types. The evolution of dust and gas disc masses inour simulations suggests a similar context for planet formation. Onthe other hand, Sellek et al. (2020) propose that photoevaporationof light dust grains can lead to a more efficient radial drift on largegrains, meaning that the dust-to-gas ratio is not greatly increased.While our current model improves previous implementations of thisratio, further modeling of the dust evolution is needed to get a betterpicture of how it affects the material available for planet formation.Another mechanism which affects our results is the long-livedstar formation process. In our previous work all the stars and discswere formed at the start of the simulations and evolved for thesame amount of time. In our current simulations, some discs aredestroyed while new stars with discs are forming constantly. By MNRAS , 1–12 (2020) Concha-Ramírez et al. the end of the simulations, the stars span ages from ∼ ∼
12 Myr, if weconsider the first stars formed in our longest-running realization asthe eldest ones. In Concha-Ramírez et al. (2019b) we show thataround 60% of discs are destroyed by photoevaporation within thefirst 100.000 years of evolution in a region of stellar density ∼
100 stars pc − . Similar rapid effects of photoevaporation are foundin Concha-Ramírez et al. (2021). While discs are also destroyedquickly in the present simulations, the fact that discs are constantlybeing added to the region helps keep the disc number and discmasses higher. Massive, radiating stars also do not form all at thesame time, and because of the way we implement the star formationprocess (see Section 2.1), they tend to form later in our simulationsthan low mass stars. This also gives discs more time to evolve, andthe older discs will be immersed in lower radiation fields duringmost of their life than the younger discs.The gas present in the clusters during the star formation processalso affects the final masses of the discs in our simulations. Whilewe do not take into account the dampening effects of the gas on radi-ation, the gas in our simulations does have important consequencesfor the dynamics of the regions. As can be seen in Figure 3, after starformation ends and the gas is expelled, the clusters quickly expandto regain virial equilibrium. This expansion leads to a decrease inthe local stellar densities and an increase in the distance betweenstars, both of which reduce the amount of radiation received by thediscs. Considering the absorption of radiation by the gas resultsin a protective effect for the discs, with the gas actively shieldingthem from external photoevaporation. The presence of gas is animportant point determined by Winter et al. (2019a) to reproducethe disc population of Cygnus OB2, and to explain the so-called‘proplyd lifetime problem’ discussed below. Observational surveysalso suggest the presence of gas being protective for the discs (e.g.van Terwisga et al. 2020). Our model shows that the effects of thegas over the dynamics of the regions could also have an importanteffect in the final disc fractions and disc mass distributions.Our present results suggest that the strong trends seen inConcha-Ramírez et al. (2019b) and Concha-Ramírez et al. (2021)of disc masses diminishing in time due to external photoevaporationmight be a temporary effect. The low mass disc distributions canbe “washed-out” in time if the star formation process lasts for anextended period of time. Massive discs that are born late in the starformation process, that is, shortly before the gas in our simulationsis expelled, will be protected by the effects of external photoe-vaporation by the expansion of the clusters as they regain virialequilibrium. In our simulations, the moment when a disc is born iscrucial for determining its long term survival and planet-formingpotential, even more so than its initial mass.The ‘proplyd lifetime problem’ is the name given to the fact thatcircumstellar discs are observed in regions were external radiationshould have evaporated them already, in particular in the ONC. Discsobserved within 0 . 𝜃 𝐶 Ori are exposed tostrong radiation fields ( ∼ G , e.g. O’dell & Wen 1994). The factthat there are still discs observed in the vicinity of the star suggestthat either the discs were initially extremely massive ( (cid:38) (cid:12) ), orthat 𝜃 𝐶 Ori is younger than 0 . ∼ We perform simulations of molecular cloud collapse and star for-mation. We begin with molecular clouds of mass 10 M (cid:12) and astar formation efficiency of 30%. The star formation process endswhen the star formation efficiency has been reached. The leftovergas is then removed instantaneously and the simulations continue.Stars with masses M ∗ ≤ . (cid:12) have circumstellar discs, and mas-sive stars (M ∗ > . (cid:12) ) emit UV radiation. The discs are subjectto internal photoevaporation, external photoevaporation, dynamicaltruncations, and viscous evolution. We also implement a model forthe evolution of dust inside the discs. We run the simulations for2 Myr after the last star has formed. We conclude that:1. Extended periods of star formation allow for relatively mas-sive (M gas ∼
100 M
𝐽𝑢 𝑝 , M dust ∼
100 M ⊕ ) discs to survivethe effects of photoevaporation for at least 2 Myr.2. The discs that make it to the end of our simulations are theones that were born later (after ∼ dust / M gas ratio increases with time. This might allow discsto keep enough mass to form rocky planets, even when de-pleted of gas.4. The presence of gas plays an important role in the survival ofcircumstellar discs, not just because it can absorb radiationfrom massive stars but also because it affects the dynamics ofthe regions.5. Extended periods of star formation can allow for massivediscs to be detected in areas with strong radiation fields.Our results show that, while photoevaporation is an important MNRAS000
100 M ⊕ ) discs to survivethe effects of photoevaporation for at least 2 Myr.2. The discs that make it to the end of our simulations are theones that were born later (after ∼ dust / M gas ratio increases with time. This might allow discsto keep enough mass to form rocky planets, even when de-pleted of gas.4. The presence of gas plays an important role in the survival ofcircumstellar discs, not just because it can absorb radiationfrom massive stars but also because it affects the dynamics ofthe regions.5. Extended periods of star formation can allow for massivediscs to be detected in areas with strong radiation fields.Our results show that, while photoevaporation is an important MNRAS000 , 1–12 (2020) ircumstellar discs in young star-forming regions process for the depletion of disc masses, other factors such as themorphology of the star-forming regions and the length of the star for-mation process can allow for circumstellar discs to survive for longperiods of time, and to remain massive enough to form planets. Thiscould help explain the great variety of disc mass distributions ob-served in star-forming regions of different ages and configurations,and could also be a factor in solving the ‘proplyd-lifetime problem’,where massive discs are observed in regions of high stellar densityand background radiation. DATA AVAILABILITY
The data underlying this article are available at https://doi.org/10.5281/zenodo.4436316
SOFTWARE
The present works makes use of the following software: AMUSE(Portegies Zwart et al. 2013; Pelupessy et al. 2013), Fi (Pelupessyet al. 2004), ph4, Bridge (Fujii et al. 2007; Portegies Zwart et al.2020), VADER (Krumholz & Forbes 2015), numpy (Van Der Waltet al. 2011), scipy (Virtanen et al. 2019), matplotlib (Hunter2007), and makecite (Price-Whelan et al. 2018).
ENERGY CONSUMPTION
The simulations presented in this work took ∼
10 days to run andused 10 cores each. Accounting for only the 6 runs shown on thispaper (not considering the many test runs), this amounts to 14.400CPU hours. Considering 12 Wh for a CPU, this results in ∼ ∼
611 kg of CO . Estimating the number oftests to be of the same order of magnitude of the final runs, the totalCO output of the computations performed for this paper are ∼ REFERENCES
Adams F. C., 2010, ARA&A, 48, 47Adams F. C., Hollenbach D., Laughlin G., Gorti U., 2004, ApJ, 611, 360Adams F. C., Proszkow E. M., Fatuzzo M., Myers P. C., 2006, ApJ, 641, 504Anderson K. R., Adams F. C., Calvet N., 2013, ApJ, 774, 9Ansdell M., et al., 2016, ApJ, 828, 46Ansdell M., Williams J. P., Manara C. F., Miotello A., Facchini S., van derMarel N., Testi L., van Dishoeck E. F., 2017, AJ, 153, 240Ansdell M., et al., 2020, arXiv:2010.00012 [astro-ph]Armitage P. J., 2000, A&A, 362, 968Ballesteros-Paredes J., et al., 2020, Space Science Reviews, 216, 76Balog Z., Muzerolle J., Rieke G. H., Su K. Y. L., Young E. T., Megeath S. T.,2007, ApJ, 660, 1532Bate M. R., Bonnell I. A., Bromm V., 2003, Monthly Notices of the RoyalAstronomical Society, 339, 577Birnstiel T., Dullemond C. P., Brauer F., 2010, A&A, 513, A79Bohlin R. C., Savage B. D., Drake J. F., 1978, ApJ, 224, 132Bontemps S., et al., 2001, ] 10.1051/0004-6361:20010474, 372, 173Breslau A., Steinhausen M., Vincke K., Pfalzner S., 2014, A&A, 565, A130Cabrit S., Pety J., Pesenti N., Dougados C., 2006, A&A, 452, 897Carpenter J. M., Mamajek E. E., Hillenbrand L. A., Meyer M. R., 2006,ApJ, 651, L49Cartwright A., Whitworth A. P., 2004, MNRAS, 348, 589Casertano S., Hut P., 1985, ApJ, 298, 80 Chen C.-Y., Li Z.-Y., King P. K., Fissel L. M., 2017, The AstrophysicalJournal, 847, 140Clarke C. J., 2007, MNRAS, 376, 1350Clarke C. J., Pringle J. E., 1991, Monthly Notices of the Royal AstronomicalSociety, 249, 588Clarke C. J., Pringle J. E., 1993, Monthly Notices of the Royal AstronomicalSociety, 261, 190Concha-Ramírez F., Vaher E., Portegies Zwart S., 2019a, MNRAS, 482,732Concha-Ramírez F., Wilhelm M. J. C., Portegies Zwart S., Haworth T. J.,2019b, MNRAS, 490, 5678Concha-Ramírez F., Wilhelm M. J. C., Zwart S. P., van Terwisga S. E.,Hacar A., 2021, Monthly Notices of the Royal Astronomical Society,501, 1782Cuello N., et al., 2018, arXiv:1812.00961 [astro-ph]Da Rio N., Robberto M., Soderblom D. R., Panagia N., Hillenbrand L. A.,Palla F., Stassun K. G., 2010, The Astrophysical Journal, 722, 1092Da Rio N., Robberto M., Hillenbrand L. A., Henning T., Stassun K. G.,2012, ApJ, 748, 14Eisner J. A., Carpenter J. M., 2006, ApJ, 641, 1162Eisner J. A., Plambeck R. L., Carpenter J. M., Corder S. A., Qi C., WilnerD., 2008, The Astrophysical Journal, 683, 304Eisner J. A., et al., 2018, ApJ, 860, 77Elmegreen B. G., Efremov Y., Pudritz R. E., Zinnecker H., 2000, Protostarsand Planets IV, p. 179Facchini S., Clarke C. J., Bisbas T. G., 2016, MNRAS, 457, 3593Falgarone E., Phillips T. G., 1991, 147, 119Falgarone E., Phillips T. G., Walker C. K., 1991, The Astrophysical Journal,378, 186Fang M., et al., 2012, A&A, 539, A119Fatuzzo M., Adams F. C., 2008, ApJ, 675, 1361Flaccomio E., Micela G., Sciortino S., 2012, Astronomy and Astrophysics,548, A85Font A. S., McCarthy I. G., Johnstone D., Ballantyne D. R., 2004, ApJ, 607,890Fujii M., Iwasawa M., Funato Y., Makino J., 2007, Publications of theAstronomical Society of Japan, 59, 1095Galli P. A. B., Bouy H., Olivares J., Miret-Roig N., Sarro L. M., BarradoD., Berihuete A., Brandner W., 2020, Astronomy & Astrophysics, 634,A98Gorti U., Hollenbach D., 2009, ApJ, 690, 1539Gorti U., Dullemond C. P., Hollenbach D., 2009, ApJ, 705, 1237–1251Guarcello M. G., et al., 2016, arXiv:1605.01773 [astro-ph]Habing H. J., 1968, Bulletin of the Astronomical Institutes of the Nether-lands, 19, 421Hacar A., Tafalla M., Forbrich J., Alves J., Meingast S., Grossschedl J.,Teixeira P. S., 2018, A&A, 610, A77Hansen C. E., Klein R. I., McKee C. F., Fisher R. T., 2012, The AstrophysicalJournal, 747, 22Hartmann L., 2002, The Astrophysical Journal, 578, 914Haworth T. J., Boubert D., Facchini S., Bisbas T. G., Clarke C. J., 2016,MNRAS, 463, 3616Haworth T. J., Facchini S., Clarke C. J., Cleeves L. I., 2017, MNRAS, 468,L108Haworth T. J., Facchini S., Clarke C. J., Mohanty S., 2018a, MNRAS, 475,5460Haworth T. J., Clarke C. J., Rahman W., Winter A. J., Facchini S., 2018b,MNRAS, 481, 452Heyer M., Goldsmith P. F., Yıldız U. A., Snell R. L., Falgarone E., PinedaJ. L., 2016, Monthly Notices of the Royal Astronomical Society, 461,3918Hillenbrand L. A., Hartmann L. W., 1998, ApJ, 492, 540Hollenbach D. J., Yorke H. W., Johnstone D., 2000, Protostars and PlanetsIV, p. 401Hunter J. D., 2007, Computing in Science & Engineering, 9, 90Jílková L., Portegies Zwart S., Pijloo T., Hammer M., 2015, MNRAS, 453,3157MNRAS , 1–12 (2020) Concha-Ramírez et al.
Jílková L., Hamers A. S., Hammer M., Portegies Zwart S., 2016, MNRAS,457, 4218Johnstone D., Hollenbach D., Bally J., 1998, ApJ, 499, 758Kim J. S., Clarke C. J., Fang M., Facchini S., 2016, ApJ, 826, L15Kraus A. L., Hillenbrand L. A., 2008, The Astrophysical Journal Letters,686, L111Kroupa P., 2001, MNRAS, 322, 231Krumholz M. R., Forbes J. C., 2015, Astronomy and Computing, 11, 1Lada C. J., Lada E. A., 2003, ARA&A, 41, 57Larson R. B., 1995, Monthly Notices of the Royal Astronomical Society,272, 213Luhman K. L., 2007, ApJS, 173, 104Luhman K. L., 2018, The Astronomical Journal, 156, 271Luhman K. L., Esplin T. L., 2020, The Astronomical Journal, 160, 44Luhman K. L., Mamajek E. E., 2012, ApJ, 758, 31Luhman K. L., Esplin T. L., Loutrel N. P., 2016, ApJ, 827, 52Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603Mac Low M.-M., Klessen R. S., 2004, Reviews of Modern Physics, 76, 125Manara C. F., et al., 2020, arXiv:2004.14232 [astro-ph]Mann R. K., Williams J. P., 2010, ApJ, 725, 430Mann R. K., et al., 2014, ApJ, 784, 82Neuhäuser R., Forbrich J., 2008, Handbook of Star Forming Regions, Vol-ume II, p. 735O’dell C. R., 1998, AJ, 115, 263O’dell C. R., Wen Z., 1994, MNRAS, 436, 194Offner S. S. R., Arce H. G., 2015, The Astrophysical Journal, 811, 146Owen J. E., Ercolano B., Clarke C. J., Alexander R. D., 2010, MNRAS, 401,1415Owen J. E., Clarke C. J., Ercolano B., 2012, MNRAS, 422, 1880Parker R. J., 2014, Monthly Notices of the Royal Astronomical Society, 445,4037Parker R. J., Alves de Oliveira C., 2017, Monthly Notices of the RoyalAstronomical Society, 468, 4340Parker R. J., Meyer M. R., 2012, MNRAS, 427, 637Pascucci I., et al., 2016, ApJ, 831, 125Pelupessy F. I., van der Werf P. P., Icke V., 2004, Astronomy and Astro-physics, 422, 55Pelupessy F. I., van Elteren A., de Vries N., McMillan S. L. W., Drost N.,Portegies Zwart S. F., 2013, A&A, 557, A84Petigura E. A., Howard A. W., Marcy G. W., 2013, Proceedings of theNational Academy of Science, 110, 19273Pfalzner S., 2003, The Astrophysical Journal, 592, 986Pfalzner S., Vogel P., Scharwächter J., Olczak C., 2005a, Astronomy andAstrophysics, 437, 967Pfalzner S., Umbreit S., Henning T., 2005b, The Astrophysical Journal, 629,526Pfalzner S., Tackenberg J., Steinhausen M., 2008, A&A, 487, L45Picogna G., Ercolano B., Owen J. E., Weber M. L., 2019, MNRAS, 487,691Portegies Zwart S. F., 2016, MNRAS, 457, 313Portegies Zwart S., 2020, Nature Astronomy, 4, 819Portegies Zwart S., McMillan S. L. W., van Elteren E., Pelupessy I., de VriesN., 2013, Computer Physics Communications, 183, 456Portegies Zwart S., Pelupessy I., Mart\’inez-Barbosa C., van Elteren A.,McMillan S., 2020, Communications in Nonlinear Science and Numer-ical Simulations, 85, 105240Price-Whelan Mechev A., jumeroag 2018, adrn/makecite: v0.2, Zenodo,doi:10.5281/zenodo.1343299Punzo D., Capuzzo-Dolcetta R., Portegies Zwart S., 2014, Monthly Noticesof the Royal Astronomical Society, 444, 2808Reche R., Beust H., Augereau J.-C., 2009, A&A, 493, 661Rodriguez-Merino L. H., Chavez M., Bertone E., Buzzoni A., 2005, ApJ,626, 411Rodriguez J. E., et al., 2018, ApJ, 859, 150Rosotti G. P., Dale J. E., de Juan Ovelar M., Hubber D. A., Kruijssen J.M. D., Ercolano B., Walch S., 2014, MNRAS, 441, 2094Rosotti G. P., Tazzari M., Booth R. A., Testi L., Lodato G., Clarke C., 2019,MNRAS, 486, 4829 Sacco G. G., et al., 2017, Astronomy and Astrophysics, 601, A97Scally A., Clarke C., 2001, MNRAS, 325, 449Scalo J., 1990, ] 10.1007/978-94-009-0605-1-12, 162, 151Sellek A. D., Booth R. A., Clarke C. J., 2020, Monthly Notices of the RoyalAstronomical Society, 492, 1279Simon M., 1997, The Astrophysical Journal Letters, 482, L81Störzer H., Hollenbach D., 1999, ApJ, 515, 669Tobin J. J., et al., 2020, ApJ, 890, 130Trapman L., Rosotti G., Bosman A. D., Hogerheijde M. R., van DishoeckE. F., 2020, arXiv:2005.11330 [astro-ph]Tritsis A., Tassis K., 2016, Monthly Notices of the Royal AstronomicalSociety, 462, 3602Van Der Walt S., Colbert S. C., Varoquaux G., 2011, Computing in Science& Engineering, 13, 22Vicente S. M., Alves J., 2005, A&A, 441, 195Vincke K., Pfalzner S., 2016, ApJ, 828, 48Vincke K., Pfalzner S., 2018, ApJ, 868, 1Vincke K., Breslau A., Pfalzner S., 2015, A&A, 577, A115Virtanen P., et al., 2019, arXiv:1907.10121 [physics]Wall J. E., McMillan S. L. W., Mac Low M.-M., Klessen R. S., Porte-gies Zwart S., 2019, arXiv e-prints, p. arXiv:1901.01132Williams J. P., Best W. M. J., 2014, ApJ, 788, 59Williams J. P., Cieza L. A., 2011, ARA&A, 49, 67Williams J. P., Blitz L., McKee C. F., 2000, Protostars and Planets IV, p. 97Winter A. J., Clarke C. J., Rosotti G., Booth R. A., 2018a, MNRAS, 475,2314Winter A. J., Clarke C. J., Rosotti G., Ih J., Facchini S., Haworth T. J., 2018b,MNRAS, 478, 2700Winter A. J., Booth R. A., Clarke C. J., 2018c, MNRAS, 479, 5522Winter A. J., Clarke C. J., Rosotti G. P., 2019a, MNRAS, 485, 1489Winter A. J., Clarke C. J., Rosotti G. P., Hacar A., Alexander R., 2019b,MNRAS, 490, 5478Winter A. J., Kruijssen J. M. D., Chevance M., Keller B. W., LongmoreS. N., 2020, MNRAS, 491, 903Wright N. J., Drake J. J., Drew J. E., Vink J. S., 2010, ApJ, 713, 871Wright N. J., Parker R. J., Goodwin S. P., Drake J. J., 2014, Monthly Noticesof the Royal Astronomical Society, 438, 639van Terwisga S. E., Hacar A., van Dishoeck E. F., 2019, A&A, 628, A85van Terwisga S. E., et al., 2020, arXiv:2004.13551 [astro-ph]This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000