https://orcid.org/0000-0002-0944-0307
ABSTRACT
The evolution of circumstellar discs is influenced by their surroundings. The relevant processes include external photoevaporation due to nearby stars and dynamical truncations. The impact of these processes on disc populations depends on the star-formation history and on the dynamical evolution of the region. Since star-formation history and the phase-space characteristics of the stars are important for the evolution of the discs, we start simulating the evolution of the star cluster with the results of molecular cloud collapse simulations. In the simulation, we form stars with circumstellar discs, which can be affected by different processes. Our models account for the viscous evolution of the discs, internal and external photoevaporation of gas, external photoevaporation of dust, and dynamical truncations. All these processes are resolved together with the dynamical evolution of the cluster, and the evolution of the stars. An extended period of star formation, lasting for at least 2 Myr, results in some discs being formed late. These late-formed discs have a better chance of survival because the cluster gradually expands with time, and a lower local stellar density reduces the effects of photoevaporation and dynamical truncation. Late formed discs can then be present in regions of high UV radiation, solving the proplyd lifetime problem. We also find a considerable fraction of discs that lose their gas content but remain sufficiently rich in solids to be able to form a rocky planetary system.
1 INTRODUCTION
Circumstellar discs are a natural consequence of the star formation process and emerge within the first 104 yr after star formation (Williams & Cieza 2011). The star formation environment, rich in gas and newly formed stars, can affect the evolution of the discs. The imprint of this environment on the young discs affects their potential to form planets and the topology of the eventual planetary system.
There are several ways in which the environment can influence the evolution of circumstellar discs. Close encounters between circumstellar discs and stellar fly-bys can affect the discs’ size, mass, and surface density. Close encounters can remove mass from the outskirts of the discs, decreasing both their mass and radius (e.g. Clarke & Pringle 1991, 1993; Pfalzner, Umbreit & Henning 2005b; Breslau et al. 2014; Vincke, Breslau & Pfalzner 2015; Portegies Zwart 2016; Vincke & Pfalzner 2016; Cuello et al. 2018; Vincke & Pfalzner 2018; Winter et al. 2018a; Concha-Ramírez, Vaher & Portegies Zwart 2019a). Several numerical implementations of this process have shown that close encounters can lead to a less steep disc’s surface density (Rosotti et al. 2014), the formation of spiral arms and other structures (Pfalzner 2003; Pfalzner et al. 2005a), accretion bursts on to the host star (Pfalzner, Tackenberg & Steinhausen 2008), and exchange of mass between discs (Pfalzner et al. 2005b; Jílková et al. 2016). Observational evidence for the effects of 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 A and argue that a recent fly-by might have caused it. Reche, Beust & Augereau (2009) demonstrate that the spiral arms observed in the disc of the triple star system HD 141569 might be the result of a fly-by.
The observed structure in GW Ori’s circum-triple disc (Smallwood et al. 2021) may have been caused either by resonances in dynamics or by planet formation. Observations by Rodriguez et al. (2018) reveal newly detected tidal streams in RW Aur A, and they propose that these might be the result of many subsequent close encounters. Winter, Booth & Clarke (2018c) simulate the disc around DO Tau, which presents a tidal tail, and argue that this shape could have been caused by a close encounter with the nearby triple system HV Tau. There is also evidence that the young disc of the Solar system was affected by such an encounter. The sharp edge of the Solar system at ∼50 au could be a sign that a passing star truncated its early disc (Breslau et al. 2014; Punzo, Capuzzo-Dolcetta & Portegies Zwart 2014). The highly eccentric and inclined orbits of the Sednitos, a group of 13 detected planetoids in the outskirts of the Solar system, suggest they might have been captured from the disc of a nearby passing star (Jílková et al. 2015).
Another mechanism that can alter the evolution of circumstellar discs is photoevaporation, through which high-energy photons heat the discs’ surfaces, causing it to evaporate. The source of these photons can be the host star (internal photoevaporation) or bright stars in the vicinity (external photoevaporation). Photoevaporation is driven by far-ultraviolet (FUV), extreme ultraviolet (EUV), and X-ray photons (Johnstone, Hollenbach & Bally 1998; Adams et al. 2004). The effects of internal and external photoevaporation on circumstellar discs are distinct. Internal photoevaporation can clear areas of the disc at specific disc radii, causing the opening of gaps (Hollenbach, Yorke & Johnstone 2000; Font et al. 2004; Fatuzzo & Adams 2008; Gorti & Hollenbach 2009; Gorti, Dullemond & Hollenbach 2009; Owen et al. 2010). External photoevaporation can remove mass from all over the disc surface, but the outer regions of the discs are more vulnerable because the material is less strongly bound to the host star (Johnstone et al. 1998; Adams et al. 2004; Haworth & Clarke 2019).
Observational evidence of external photoevaporation was first obtained through the imaging of evaporating discs in the Orion nebula (O’dell & Wen 1994; O’dell 1998). These objects, now known as ‘proplyds’, are circumstellar discs immersed in the radiation fields of nearby stars. Their cometary tail-like structure reveals the ongoing mass-loss. Subsequent observations of the region showed that circumstellar disc masses decrease when close to massive stars. This effect has been observed in several regions such as the Trapezium cluster (e.g. Vicente & Alves 2005; Eisner & Carpenter 2006; Mann et 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 (Fang et al. 2012), NGC 2024 (van Terwisga et al. 2020), σ Orionis (Ansdell et al. 2017), and λ Orionis (Ansdell et al. 2020). Younger and low-mass star-forming regions such as Lupus, Taurus, Ophiuchus, and the Orion Molecular Cloud 2 tend to have higher average disc masses than denser regions such as the Orion Nebula Cluster (Eisner et al. 2008; Ansdell et al. 2016; Eisner et al. 2018; van Terwisga, Hacar & van Dishoeck 2019). van Terwisga et al. (2020) present the discovery of two distinct disc populations, in terms of mass, in the NGC 2024 region. The discs to the east of the region are embedded in a dense molecular ridge and are more massive than the discs outside the ridge, which are also closer to two OB-type stars. They propose that the difference in masses is caused by the eastern population being protected from the radiation of the nearby massive star IRS 1.
Several models have demonstrated that external photoevaporation is efficient in depleting disc masses on time-scales much shorter than their estimated lifetimes of ∼10 Myr (e.g. Scally & Clarke 2001; Adams et al. 2006; Fatuzzo & Adams 2008; Haworth et al. 2016), even in low radiation fields (Facchini, Clarke & Bisbas 2016; Kim et al. 2016; Haworth et al. 2017). Because external photoevaporation is caused by massive stars in the vicinity, the extent of its effect depends on the density of the stellar region and the number of massive stars in the surroundings. In high-density regions (N* ≳ 104 pc−3) where the stellar population follows a well-sampled IMF, the disc mass-loss rates caused by external photoevaporation are orders of magnitude higher than those caused by dynamical truncations (Winter et al. 2018b, 2019; Concha-Ramírez et al. 2019b). Concha-Ramírez et al. (2021) show that, in regions of local stellar densities N* > 100 pc−3, external photoevaporation can evaporate up to 90 per cent of circumstellar discs within 2.0 Myr. In low density regions (∼10 M⊙ pc−3), only |$\sim 60{{\ \rm per\ cent}}$| of discs are evaporated within the same time-scale. Winter et al. (2020a) model a region comparable to the central molecular zone of the Milky Way (surface density Σ0 = 103 M⊙ pc−2) and find that external photoevaporation destroys 90 per cent of circumstellar discs within 1.0 Myr. In regions of lower density (Σ0 = 12 M⊙ pc−2) they find a mean disc dispersal time-scale of 3.0 Myr. Similar results are obtained by Nicholson et al. (2019) who find that external photoevaporation destroys 50 per cent of discs within |$1.0\, \textrm {Myr}$| in regions of density ∼100 M⊙ pc−3, and within |$2.0\, \textrm {Myr}$| in regions of density ∼10 M⊙ pc−3.
While observational and numerical evidence indicate that disc masses decrease with increasing stellar density, massive discs are still observed in high-density regions. For example in the ONC, which contains discs with a mass higher than those in the proximity of the massive star θ1 Ori C. If the discs are coeval with θ1 Ori C, they should have already been dispersed by external photoevaporation, unless they were extraordinarily massive to begin with (Mdisc ≳ 1 M⊙). Alternatively, θ1 Ori C would have to be considerably younger than the ONC average (|$\lesssim 0.1\, \textrm {Myr}$|) for these discs to have survived. This is known as the ‘proplyd lifetime problem’. Störzer & Hollenbach (1999) propose that these discs are currently passing through the region’s centre, whereas they spent most time at a larger distance, where the radiation is ineffective. Scally & Clarke (2001) model external photoevaporation on a cluster similar to the ONC and find that the necessary radial orbits proposed by Störzer & Hollenbach (1999) are not dynamically plausible in such a region. Winter et al. (2019) revisit the problem and propose a solution to describe why these discs exist. Their solution consists of a combination of factors: different eras of star formation allow for massive discs to be around stars that are younger than the average of the ONC population; stars forming in subvirial states with respect to the gas potential allows young stars to migrate to the central region of the ONC; and interstellar gas protects the discs from the radiation, allowing them to live longer than expected.
The star-formation history and primordial stellar distributions in young star-forming regions are key to understanding the environment’s effects on the disc populations. The star formation process results in regions with different morphologies, with clumps and filaments likely to be present. This structure is far from the spherical, idealized initial conditions commonly used in models of star clusters. The collapse of the giant molecular clouds (GMCs) from which stars form is affected by turbulent flows (Falgarone & Phillips 1991; Falgarone, Phillips & Walker 1991) which result in filamentary, clumpy, or fractal gas substructure in the cloud (e.g. Scalo 1990; Larson 1995; Elmegreen et al. 2000; Hacar et al. 2018). The distribution of the newly formed stars follows the local densities of the gas in the molecular clouds. Star clusters would then form in a bottom-up scenario, where subclusters emerge from clumps and filaments with high star-formation efficiency that eventually merge (see Fujii & Portegies Zwart 2016). This has been determined theoretically (e.g. Elmegreen 2008; Kruijssen 2012) and suggested by recent observations (e.g. Ward, Kruijssen & Rix 2020). Other authors (e.g. Banerjee & Kroupa 2015; Kuhn et al. 2019) have proposed a monolithic or top-down formation scenario for star clusters, i.e. one where a single, highly active star-formation episode yields a compact cluster that remains embedded in its primordial gas. Once feedback processes expel this gas, the cluster expands, and loses stars. While circumstellar discs emerge during the protostellar phase (Williams & Cieza 2011), the star formation process defines the environment in which the discs are immersed at their early stages.
It is important to take a step back in time to better understand the environmental effects on circumstellar discs and study how the star formation process influences stellar densities. This work presents a model for circumstellar discs inside young star-forming regions. We adopt a relatively simple model for the star formation process, starting from the collapse of a giant molecular cloud to obtain masses, positions, and velocities of newly formed stars. These form the input for our star cluster evolution code. During the evolution, we take into account the viscous evolution of the discs, dynamical truncations, external and internal photoevaporation of gas, and external photoevaporation of dust. We evolve the discs simultaneously with the stellar dynamics and stellar evolution.
2 MODEL
We simulate the formation of massive clusters of bound stars, including circumstellar discs. These discs are subject to viscous spreading, dynamical truncations, and photoevaporation.
We model several different astrophysical processes which operate simultaneously on a range of scales: the collapse of a molecular cloud, star formation, stellar dynamics, and viscous circumstellar discs, which are affected by dynamical truncations and photoevaporation. We bring these processes together using the Astrophysical Multipurpose Software Environment, amuse (Pelupessy et al. 2013; Portegies Zwart et al. 2013). The results presented in this work are obtained through two different simulation stages: first, we simulate the collapse of a molecular cloud, including star formation. This results is a distribution of stars in space with ages, masses, and velocities. We continue to evolve these stars in the second stage, including the evolution of the circumstellar discs. This second stage encompasses the stellar dynamics, stellar evolution, viscous evolution of the discs, and photoevaporation. All the code developed for this work is available online.
2.1 Stage 1: The collapse of the molecular cloud and the formation of stars
In the first stage, we model the star-formation process. We simulate the collapse of a molecular cloud using the smoothed particle hydrodynamics (SPH) code Fi (Pelupessy, van der Werf & Icke 2004). The simulations start with a spherical gas reservoir of mass 104 M⊙ with an initial radius of |$3\, \textrm {pc}$|. The SPH simulations use 32 000 equal-mass particles, with a mass resolution of 0.3 M⊙ per particle. The softening in the simulations is ϵ = 0.05 pc. A power-law velocity spectrum P(k) ∝ k−4 was adopted to emulate large scale turbulence (here k is the wavenumber; Bate, Bonnell & Bromm 2003). These initial conditions result in a velocity dispersion in three dimensions that scales with the observed relations by Larson (1981). The mean temperature of the gas is 23 K, and the initial free-fall timescale is |$0.42\, \textrm {Myr}$|.
We model the star formation process using sink particles, which are created from regions of the cloud where the local gas density is higher than 1M⊙/ϵ3 ≃ 8000 M⊙ pc−3. Once formed, sink particles continue to accrete gas from the molecular cloud. Each of these sink particles can form several stars.
Since we aim to preserve a power-law stellar mass function similar to the one observed in the galaxy, the stars in our simulations are formed by introducing a predetermined star-formation efficiency (SFE). We set an SFE of 0.3 (Lada & Lada 2003). In reality, the integrated SFE of molecular clouds is much lower, and values around 0.3 can be found around much smaller cores (e.g. Matzner & McKee 2000; Chevance et al. 2020; Smith et al. 2020).
We implement this SFE by keeping track of the total initial mass of all stars. When this mass exceeds 30 per cent of the initial cloud mass, the SPH code is stopped, and all remaining gas and sink particles are instantaneously removed.
The star formation process begins once sink particles have accreted enough mass to sample stars randomly from the initial mass function (IMF). We base the star formation mechanism on Wall et al. (2019). We begin by sampling a random stellar mass m from a Kroupa IMF (Kroupa 2001) of 10.000 stars, with a lower limit of 0.08 M⊙ and upper limit 150 M⊙ (Wall et al. 2019). Then, the list of sinks is checked to find one that is massive enough to form the next star in the list. The first sink to satisfy this condition will be selected. We subtract the mass m from the sink, and a new star is born. If m ≤ 1.9M⊙ the star will have a circumstellar disc (see Section 2.2.2), and we subtract the mass of the star and the initial mass of the disc from the sink. The position of the newly formed star is determined by taking the position of the sink, and we add a random offset in each spatial dimension. This offset is uniformly distributed within plus and minus the sink radius. The sink radius corresponds to the distance from which the sink accretes material. By the time star formation has started, sink radii are 8 × 103 au. The velocity of the new star is set to the velocity of its birth sink.
After a sink has formed a star, we set a delay time that must pass before creating a new star. We implement this step to prevent the instantaneous conversion of gas into stars. This is typical if only hydrodynamics and self-gravity are accounted for, but observations indicate that star formation proceeds at a lower rate (Krumholz, McKee & Bland-Hawthorn 2019). Various processes can contribute to this, such as magnetic fields and stellar feedback, which we do not simulate. The delay is implemented as an exponentially decaying time-scale of
Here tff is the free-fall time-scale of the corresponding sink, t is the current model time, and td is a free parameter we fix at |$1\, \textrm {Myr}$|. It is not clear a priori how this recipe translates to a star formation rate, but we explore this in Section 4.2.
During the molecular cloud-collapse simulation, we keep track of the newly formed stars’ mass, position, velocity, and birth time. This data will then be input for the second part of the simulations, in which the discs are evolved together with the stellar dynamics. The dynamical evolution of the stars is calculated using the 4th-order Hermite integrator ph4, which is integrated together with the SPH 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). Stellar evolution was accounted for using the SeBa stellar and binary evolution package (Portegies Zwart & Verbunt 1996; Toonen et al. 2020).
2.2 Stage 2: The dynamics of stars and the evolution of their circumstellar discs
The second simulation stage begins when the first star has formed. For each star, we evolve its disc and calculate its external photoevaporation mass-loss rate as explained in Sections 2.2.1 and 2.2.2, respectively. The star formation process ends when 30 per cent of the initial cloud mass has been converted into stars. Then, the SPH code and Bridge coupling between gas and stars are stopped, in effect instantaneously removing the excess gas. From then on we only deal with the stars’ dynamics and evolution, and the processes as explained in the following sections. This second stage of the simulations run for |$2\, \textrm {Myr}$| after the last star has formed. Below we describe how each of these process is incorporated.
2.2.1 Circumstellar discs
We model circumstellar discs using the Viscous Accretion Disk [sic] Evolution Resource (VADER; Krumholz & Forbes 2015). VADER models the viscous transport of mass and angular momentum in thin, axisymmetric discs. Each disc is defined with a grid of 100 logarithmically spaced cells, spanning 0.05 to |$2000\, \textrm {au}$|. The larger outer limit allows the discs to expand without reaching the grid boundaries. The mass flow through the outer boundary of the grid is set to zero to maintain the density required to measure the disc radius. The mass flow from the disc’s inner boundary is considered accreted on to the host star. Each disc has a Keplerian rotation profile and turbulence parameter α = 5 × 10−3, which results in a viscous time-scale of |$\sim 0.1\, \textrm {Myr}$|.
We use the standard disc profile of Lynden-Bell & Pringle (1974) to establish the initial column density of the discs as:
Here rd is the initial disc radius, md is the initial disc mass. We consider the characteristic radius to be rc ≈ rd (Anderson, Adams & Calvet 2013).
For the external photoevaporation process, we keep track of the outer edge of the disc. We define the disc radius rd as the radius which encloses 99 per cent of the disc mass (Anderson et al. 2013), and set the column density outside rd to a negligible value |$\Sigma _{\mathrm{edge}} = 10^{-12}\, \mathrm{g\, cm}^{-3}$|. The mass-loss due to external photoevaporation (Section 2.2.2), as well as dynamical truncations (Section 2.2.5), causes the disc to develop a steep density profile at the outer edge. The location of the edge is insensitive to the value of Σedge, given that it is sufficiently low (Clarke 2007).
2.2.2 External photoevaporation
We calculate the mass-loss due to external photoevaporation using the Far-ultraviolet Radiation Induced Evaporation of Discs (FRIED) grid (Haworth et al. 2018b). This grid provides a pre-calculated set of mass-loss rates for discs immersed in UV radiation fields of varying strengths, from 10 G0 to 104 G0, where G0 is the FUV field measured in Habing units, 1.6 × 10−3 erg s−1 cm−3 (Habing 1968). The grid spans discs of mass ∼10−4 MJup to 102 MJup, radius from |$1\, \textrm {au}$| to |$400\, \textrm {au}$|, and host star mass from 0.05 M⊙ to 1.9 M⊙.
To stay within the boundaries of the FRIED grid, we consider only stars of M* ≤ 1.9 M⊙ to have a circumstellar disc. More massive stars are considered to generate too much radiation to have an appreciable disc. This mass distinction is for external photoevaporation calculations only; for the stellar dynamics evolution, there is no separation between these two stellar groups.
Far-ultraviolet (FUV) photons dominate in the external photoevaporation process (Armitage 2000; Adams et al. 2004; Gorti & Hollenbach 2009). We calculate the FUV radiation from the massive stars by pre-computing a relation between stellar mass and FUV luminosity using the UVBLUE spectral library (Rodriguez-Merino et al. 2005). The obtained fit is presented in fig. 2 of Concha-Ramírez et al. (2019b). We use that fit to determine the FUV radiation emitted by each massive star at every simulation time-step.
The external photoevaporation process is applied at every time-step. We first calculate the distance from every disc to every star of mass M* > 1.9 M⊙ and determine the total radiation received by each disc. We do not consider extinction due to interstellar material. We interpolate from the FRIED grid using the calculated total radiation and the disc parameters to find a photoevaporation-driven mass-loss rate for each disc. Assuming the mass-loss rate |$\dot{\mathrm{M}}$| to be constant during the current time-step, we use it to calculate the total mass-loss. This mass is subsequently removed from the disc’s outer region, removing mass from each subsequent cell until a finite amount of mass remains. External photoevaporation then results in a decrease of disc mass and disc radius.
External photoevaporation can be dominated by extreme ultraviolet (EUV) photons. This happens when a disc is closer to a radiating star than (Johnstone et al. 1998)
Here fr is the fraction of EUV photons absorbed in the ionizing flow, Φ49 = (Φi/1049)s−1 is the EUV luminosity of the source, ε is a dimensionless normalizing parameter, (ε2/(frΦ49))1/2 ≈ 4, and |$r_{d_{14}} ={r_d}/({10^{14} \mathrm{cm}})$| with rd the disc radius.
When the distance d between a disc and a massive star is d < dmin, the disc enters the EUV-dominated photoevaporation regime. The mass-loss in this case is calculated as:
with x ≈ 1.5 and ε ≈ 3 (Johnstone et al. 1998).
A disc is considered dispersed when it has lost 99 per cent of its initial gas mass (Ansdell et al. 2016) or when its mean column density drops below |$1\, \mathrm{g\, cm}^{-2}$| (Pascucci et al. 2016). After a disc is dispersed, its host star continues to be integrated by the stellar dynamics and stellar evolution codes.
2.2.3 Internal photoevaporation
Internal photoevaporation is driven by X-ray radiation (Owen et al. 2010; Owen, Clarke & Ercolano 2012). We calculate the X-ray luminosity of the stars with discs using the fit obtained by Flaccomio, Micela & Sciortino (2012) for T-Tauri stars as a function of stellar mass:
Picogna et al. (2019) calculate X-ray mass-loss profiles and mass-loss rates for a star of mass 0.7 M⊙. Owen et al. (2012) developed scaling relations that allow us to calculate these values for stars of M* ≤ 1.5M⊙. We combine the results of Picogna et al. (2019) with the scaling relations of Owen et al. (2012) to span a larger range of stellar masses. The internal photoevaporation mass-loss rate is then given by
where
is the X-ray mass-loss rate derived by Picogna et al. (2019).
The mass-loss profile takes the form
where
Here a = −0.5885, b = 4.3130, c = −12.1214, d = 16.3587, e = −11.4721, f = 5.7248, and g = −2.8562 (Picogna et al. 2019), and
is the scaling from Owen et al. (2012).
The internal photoevaporation process removes mass from the disc following the profile defined in equation (8). If a grid cell contains less mass than is prescribed to be removed, this excess is removed in the nearest outer cell. As the cells are traversed inside-out, this takes the form of inside-out disc clearing.
2.2.4 Disc dust photoevaporation
Circumstellar discs are composed of gas and dust. Initially, the gas:dust ratio is 100:1. This value is derived from the consideration that the ratio is inherited from the interstellar medium (Bohlin, Savage & Drake 1978). Grain growth might result in much lower gas:dust ratios (Williams & Best 2014), and the ratio is likely to change during the lifetime of a disc (Manara et al. 2020). Models of dust evolution and radial drift (Birnstiel, Dullemond & Brauer 2010; Rosotti et al. 2019) show that the dust-to-gas mass-ratio δ decreases with time. We introduce a simple prescription for the photoevaporation of dust inside circumstellar discs in this work.
We follow the prescription of Haworth et al. (2018a) to calculate the mass-loss rate of dust entrapped in the photoevaporation wind. This mass-loss rate is described as:
Here the initial dust-to-gas mass-ratio δ = 10−2, |$\dot{M}_{\rm gas}$| is the gas mass-loss rate (determined as explained in Sections 2.2.2 and 2.2.3), |${\nu }_{\rm th} = \sqrt{8 k_{b} T / (\pi \mu m_{H})}$| is the mean thermal speed of the gas particles, F is the solid angle subtended by the disc at the outer edge, ρg is the grain mass density (|$1 \, \mathrm{g\, cm}^{-3}$|; Facchini et al. 2016), and amin is the minimum grain size at the disc radius Rd. We assume amin = 0.01 |$\mu$|m (Facchini et al. 2016; Haworth et al. 2018a).
This model takes into account, the fraction of the dust entrained in the photoevaporation wind, and how it decreases over time due to dust growth. Mass is removed from a single scalar reservoir. The radial structure is implicitly assumed to follow the gas structure, multiplied by the dust-to-gas ratio δ ∼ 0.01, and does not account for the dust fraction enhancement due to the evaporation-resistant dust population.
This reservoir can also be depleted by other processes such as inward radial drift and the formation of pebbles, planetesimals, and planets. For this reason, we consider the dust reservoir in this model to be an upper limit for the total mass in solids with an unknown size distribution.
2.2.5 Dynamical truncations
Circumstellar discs can be truncated in encounters with other stars in the cluster. We calculate a semi-analytical truncation radius based on Adams (2010), who propose that the new radius of a disc after a truncating encounter is R′ ≈ b/3, where b is the pericentre distance of the encounter. We combine this with the mass dependence of Breslau et al. (2014) to define a truncation radius:
Here m1 and m2 are the masses of the encountering stars. We ignore the disc orientation, and the equation for truncation radius is the average truncation radius over all inclinations. We follow Portegies Zwart (2016) in defining an initial collisional radius of |$r_{\rm col} = 0.02\, \textrm {pc}$| for all stars. This value is updated to rcol = 0.5renc after every encounter, to guarantee that each encounter is only detected once within the time-step. Not all encounters result in disc truncation. If the calculated truncation radius R′ exceeds the current radius of an encountering disc, the disc is not affected by the encounter. If a disc is truncated in an encounter, we set the new radius of a disc to R′ by making the column density |$\Sigma _{\mathrm{edge}} = 10^{-12}\, \mathrm{g\, cm}^{-3}$| for every disc cell outside R′. The truncated disc then continues to expand viscously.
Dynamical encounters not only change the disc sizes and strip mass from the outskirts but can also lead to changes in the mass distribution of the discs, and mass exchange can occur between the encountering discs (Pfalzner et al. 2005b; Rosotti et al. 2014; Jílková et al. 2015; Portegies Zwart 2016). We do not consider mass exchange or changes to the mass distribution during dynamical encounters, other than truncation. When a disc is truncated in our model, all the mass outside its new radius R′ is simply lost.
2.3 Initial conditions
2.3.1 Molecular cloud
Our simulations start with a spherical cloud model of mass 104M⊙ and initial radius |$3\, \textrm {pc}$|. We use 32 000 SPH particles, which results in a resolution of 0.3 M⊙ per particle. The softening in the simulations is |$0.05\, \textrm {pc}$|. We use a power-law velocity spectrum to model large-scale turbulence (Bate et al. 2003). Each realization of a molecular cloud has a different random seed, and as a result, the substructure is different in every run. We run six realizations of the molecular cloud collapse simulations. These realizations differ only in the random seed used to determine the position and velocities of the SPH particles.
2.3.2 Circumstellar discs
We choose the initial disc radii as:
Here R′ is a constant. We choose |$R^{\prime } = 30\, \textrm {au}$|, which results in initial disc radii between |$\sim 5\, \textrm {au}$| and |$\sim 40\, \textrm {au}$|. This is in agreement with observations that suggest that young circumstellar discs are generally quite compact (radii around 20 to |$50\, \textrm {au}$|; Tobin et al. 2020; Trapman et al. 2020).
The initial mass of the discs:
This yields initial disc masses ranging from ∼8 MJup to ∼200 MJup.
3 RESULTS
3.1 Star formation and cluster evolution
In Fig. 1 we illustrate the evolution of the molecular cloud collapse and star formation process in one of our simulations. In Fig. 2 we show the number of stars in time for each simulations. In Table 1 we present the final number of stars in each run. The mean number of stars created in six runs is 5031 ± 198.
The evolution of the molecular cloud collapse and star formation process in run #4. The second panel shows the moment before the first stars form, and the centre panel shows the region during star formation. The fourth panel shows the moment just before gas expulsion, and the rightmost panel shows the region close to the end of the simulation.
Number of stars in time for each simulation run.
Run number, final number of stars (N*), time of the end of star formation (|$\mathit{ t}^{\mathrm{SF}}_{\mathrm{end}}$|), Q parameter, and fractal dimension (Fd) for our simulation results.
| Region | N* | |$\mathit{ t}^{\mathrm{SF}}_{\mathrm{end}}$| (Myr) | Qend | Fd |
|---|---|---|---|---|
| Run #1 | 5198 | 4.20 | 1.11 | 1.8 |
| Run #2 | 5242 | 4.75 | 0.61 | 1.0 |
| Run #3 | 5071 | 3.05 | 1.06 | 1.8 |
| Run #4 | 4654 | 3.32 | 0.96 | 1.6 |
| Run #5 | 4914 | 3.31 | 1.06 | 1.6 |
| Run #6 | 5106 | 5.51 | 1.07 | 1.7 |
| Region | N* | |$\mathit{ t}^{\mathrm{SF}}_{\mathrm{end}}$| (Myr) | Qend | Fd |
|---|---|---|---|---|
| Run #1 | 5198 | 4.20 | 1.11 | 1.8 |
| Run #2 | 5242 | 4.75 | 0.61 | 1.0 |
| Run #3 | 5071 | 3.05 | 1.06 | 1.8 |
| Run #4 | 4654 | 3.32 | 0.96 | 1.6 |
| Run #5 | 4914 | 3.31 | 1.06 | 1.6 |
| Run #6 | 5106 | 5.51 | 1.07 | 1.7 |
Run number, final number of stars (N*), time of the end of star formation (|$\mathit{ t}^{\mathrm{SF}}_{\mathrm{end}}$|), Q parameter, and fractal dimension (Fd) for our simulation results.
| Region | N* | |$\mathit{ t}^{\mathrm{SF}}_{\mathrm{end}}$| (Myr) | Qend | Fd |
|---|---|---|---|---|
| Run #1 | 5198 | 4.20 | 1.11 | 1.8 |
| Run #2 | 5242 | 4.75 | 0.61 | 1.0 |
| Run #3 | 5071 | 3.05 | 1.06 | 1.8 |
| Run #4 | 4654 | 3.32 | 0.96 | 1.6 |
| Run #5 | 4914 | 3.31 | 1.06 | 1.6 |
| Run #6 | 5106 | 5.51 | 1.07 | 1.7 |
| Region | N* | |$\mathit{ t}^{\mathrm{SF}}_{\mathrm{end}}$| (Myr) | Qend | Fd |
|---|---|---|---|---|
| Run #1 | 5198 | 4.20 | 1.11 | 1.8 |
| Run #2 | 5242 | 4.75 | 0.61 | 1.0 |
| Run #3 | 5071 | 3.05 | 1.06 | 1.8 |
| Run #4 | 4654 | 3.32 | 0.96 | 1.6 |
| Run #5 | 4914 | 3.31 | 1.06 | 1.6 |
| Run #6 | 5106 | 5.51 | 1.07 | 1.7 |
In Fig. 3 we show the evolution of the cluster half-mass radius in time for each of our simulations. The solid lines show the half-mass radius while star formation is still ongoing, whereas the dotted lines follow the radius after all stars have formed and gas has been removed from the clusters. All regions initially expand gradually at a rate of 1 to 2 pc Myr−1. After gas expulsion this expansion accelerates to a rate of 3 to 10 pc Myr−1.
Cluster half-mass radius of the simulations in time. The solid lines correspond to the half-mass radius while star formation is still ongoing and the dotted lines after it has ended.
To quantify the spatial distribution of the stars resulting from the molecular cloud collapse simulations, we look at the Q parameter of the minimum spanning tree (Cartwright & Whitworth 2004) and the fractal dimension in each region at the end of star formation. The Q parameter is calculated as:
Here |$\overline{m}$| is the mean length of the minimum spanning tree, and |$\overline{s}$| is the mean separation between the stars. Regions with Q > 0.8 are smooth and centrally concentrated, while values of Q < 0.8 correspond to regions with substructure.
In Fig. 4 we show the Q parameter at time for each of our simulations, along with values for several observed regions. The Q parameter of the simulations is calculated from a 2D projection of the stellar distances, and considers only stars with masses M* > 0.5 M⊙. In Table 2 we summarize the values for Q and the estimated ages for each region, along with the corresponding references. The simulation results span a range of Q ∼ 0.5–0.8 when star formation starts and evolve to Q ∼ 0.9–1.1.
Q parameter of our simulations in time, and values for observed star-forming regions. The Q parameter in our simulations considers only stars with masses M* > 0.5 M⊙. The Q parameters are computed at times corresponding to the ages of the observed star-forming regions and the end of the simulation, and shown as points connected by linear dashed lines to guide the eye.
Region name, number of stars (N*), age, Q parameter, and fractal dimension (|$F\_d$|) for our simulation results and observed regions.
| Region | N* | Age (Myr) | Q | Fd |
|---|---|---|---|---|
| Corona Australis (CrA) | ∼313q | ∼1.0a | 0.32b, 0.38a | – |
| NGC 1333 | ∼200c | ∼1.0c | 0.89d | – |
| ONC | ∼1000e | ∼1.0e | 0.87e, 0.94a | – |
| Taurus | ∼438s | ∼1.0f | 0.47f | 1.5f, 1.02 ± 0.04g, 1.049 ± 0.007h, 1.5 ± 0.2i |
| Trapezium | ∼1000e | ∼1.0e | – | 1.5 ± 0.2i |
| Ophiuchus | 199f | 1.6 ± 1.4j | 0.85f, 0.58a | 1.5 ± 0.2i |
| Chamaeleon I | 120p | 2.5 ± 0.5k | 0.67f, 0.71a, 0.80 ± 0.08p | 2.25f |
| Cygnus OB2 | ∼2700l | 4.0 ± 1.0l | 0.45 ± 0.05m | – |
| IC 348 | ∼500c | 4.0 ± 2.0c | 0.98d | – |
| Upper Scorpio | ∼1761r | 8.0 ± 3.0n | 0.88h, 0.75a | 0.69 ± 0.09h |
| Region | N* | Age (Myr) | Q | Fd |
|---|---|---|---|---|
| Corona Australis (CrA) | ∼313q | ∼1.0a | 0.32b, 0.38a | – |
| NGC 1333 | ∼200c | ∼1.0c | 0.89d | – |
| ONC | ∼1000e | ∼1.0e | 0.87e, 0.94a | – |
| Taurus | ∼438s | ∼1.0f | 0.47f | 1.5f, 1.02 ± 0.04g, 1.049 ± 0.007h, 1.5 ± 0.2i |
| Trapezium | ∼1000e | ∼1.0e | – | 1.5 ± 0.2i |
| Ophiuchus | 199f | 1.6 ± 1.4j | 0.85f, 0.58a | 1.5 ± 0.2i |
| Chamaeleon I | 120p | 2.5 ± 0.5k | 0.67f, 0.71a, 0.80 ± 0.08p | 2.25f |
| Cygnus OB2 | ∼2700l | 4.0 ± 1.0l | 0.45 ± 0.05m | – |
| IC 348 | ∼500c | 4.0 ± 2.0c | 0.98d | – |
| Upper Scorpio | ∼1761r | 8.0 ± 3.0n | 0.88h, 0.75a | 0.69 ± 0.09h |
Notes. References: aParker (2014); bNeuhäuser & Forbrich (2008); cLuhman, Esplin & Loutrel (2016); dParker & Alves de Oliveira (2017); eHillenbrand & Hartmann (1998); fCartwright & Whitworth (2004); gHartmann (2002); hKraus & Hillenbrand (2008); iSimon (1997); jBontemps et al. (2001); kLuhman (2007); lWright et al. (2010); mWright et al. (2014); nCarpenter et al. (2006); oLuhman & Mamajek (2012); pSacco et al. (2017); qGalli et al. (2020); rLuhman & Esplin (2020); sLuhman (2018).
Region name, number of stars (N*), age, Q parameter, and fractal dimension (|$F\_d$|) for our simulation results and observed regions.
| Region | N* | Age (Myr) | Q | Fd |
|---|---|---|---|---|
| Corona Australis (CrA) | ∼313q | ∼1.0a | 0.32b, 0.38a | – |
| NGC 1333 | ∼200c | ∼1.0c | 0.89d | – |
| ONC | ∼1000e | ∼1.0e | 0.87e, 0.94a | – |
| Taurus | ∼438s | ∼1.0f | 0.47f | 1.5f, 1.02 ± 0.04g, 1.049 ± 0.007h, 1.5 ± 0.2i |
| Trapezium | ∼1000e | ∼1.0e | – | 1.5 ± 0.2i |
| Ophiuchus | 199f | 1.6 ± 1.4j | 0.85f, 0.58a | 1.5 ± 0.2i |
| Chamaeleon I | 120p | 2.5 ± 0.5k | 0.67f, 0.71a, 0.80 ± 0.08p | 2.25f |
| Cygnus OB2 | ∼2700l | 4.0 ± 1.0l | 0.45 ± 0.05m | – |
| IC 348 | ∼500c | 4.0 ± 2.0c | 0.98d | – |
| Upper Scorpio | ∼1761r | 8.0 ± 3.0n | 0.88h, 0.75a | 0.69 ± 0.09h |
| Region | N* | Age (Myr) | Q | Fd |
|---|---|---|---|---|
| Corona Australis (CrA) | ∼313q | ∼1.0a | 0.32b, 0.38a | – |
| NGC 1333 | ∼200c | ∼1.0c | 0.89d | – |
| ONC | ∼1000e | ∼1.0e | 0.87e, 0.94a | – |
| Taurus | ∼438s | ∼1.0f | 0.47f | 1.5f, 1.02 ± 0.04g, 1.049 ± 0.007h, 1.5 ± 0.2i |
| Trapezium | ∼1000e | ∼1.0e | – | 1.5 ± 0.2i |
| Ophiuchus | 199f | 1.6 ± 1.4j | 0.85f, 0.58a | 1.5 ± 0.2i |
| Chamaeleon I | 120p | 2.5 ± 0.5k | 0.67f, 0.71a, 0.80 ± 0.08p | 2.25f |
| Cygnus OB2 | ∼2700l | 4.0 ± 1.0l | 0.45 ± 0.05m | – |
| IC 348 | ∼500c | 4.0 ± 2.0c | 0.98d | – |
| Upper Scorpio | ∼1761r | 8.0 ± 3.0n | 0.88h, 0.75a | 0.69 ± 0.09h |
Notes. References: aParker (2014); bNeuhäuser & Forbrich (2008); cLuhman, Esplin & Loutrel (2016); dParker & Alves de Oliveira (2017); eHillenbrand & Hartmann (1998); fCartwright & Whitworth (2004); gHartmann (2002); hKraus & Hillenbrand (2008); iSimon (1997); jBontemps et al. (2001); kLuhman (2007); lWright et al. (2010); mWright et al. (2014); nCarpenter et al. (2006); oLuhman & Mamajek (2012); pSacco et al. (2017); qGalli et al. (2020); rLuhman & Esplin (2020); sLuhman (2018).
Run #2 forms a bit of an exception is several regards. Part of its deviating values (see Table 1), in Q but also in the fractal dimension probably result from an event that occurred some time during the growth of two massive sinks. These two sinks are ejected from the cluster at a velocity of 200 and 600 km s−1, possibly resulting from a three-body interaction. By the end of the simulation both sinks have produced ∼100 stars that co-move with their parent sinks. We did not want to remove this run from the simulations, even though we suspect that this ejection is the result of a numerical error in the hydrodynamics solver. We have to note there that the combined contributions of all the processes in our simulations make it hard to check for energy conservation, because in the combination of irradiation, internal disc evolution, stellar evolution, and gravitational dynamics such conservation laws do not apply locally. The main cluster, from which the two smaller clumps are ejected, turns out rather usual, and we decided to take this cluster into account in the further analysis.
In observations, the Q parameter can greatly vary depending on stellar membership. Regions with lower Q, such as Corona Australis (Parker 2014), Cygnus OB2 (Wright et al. 2014), and Taurus (Cartwright & Whitworth 2004) are highly substructured. In observations of star-forming regions, Q might vary depending on membership uncertainty (Parker & Meyer 2012). This leads to regions with more than one Q value, such as Corona Australis, Chamaeleon, the ONC, Ophiuchus, and Upper Scorpio.
As can be seen in Fig. 4, the final shapes of the clusters generated by our simulations are smoother than those of observed clusters. The smoothness of our simulated clusters might be caused by the absence of stellar feedback (as was also argued in Mac Low & Klessen 2004; Hansen et al. 2012; Offner & Arce 2015, see also Section 4).
Another way to quantify the structure of a star-forming region is by measuring its fractal dimension, Fd. The fractal dimension is a measurement of the clumpiness of a region. Low values of Fd (≲ 1.5) signify regions with important substructure. Higher values (∼2.0 to 3.0) mean that the regions are smoother (Goodwin & Whitworth 2004; de La Fuente Marcos & de La Fuente Marcos 2006). In Fig. 5 we show the evolution of the fractal dimension in time for each simulation. Initially they span a range from Fd ∼ 0.9–1.8, and first increase before decreasing to a general range of Fd ∼ 1.5–1.8. Again, run #2 is an exception due to its runaway clusters.
Fractal dimension of our simulations in time, and values for observed star-forming regions. The fractal dimension in our simulations considers only stars with masses M* > 0.5 M⊙. The fractal dimensions are computed at times corresponding to the ages of the observed star-forming regions and the end of the simulation, and shown as points connected by linear dashed lines to guide the eye.
3.2 Disc masses
In Figs 6 and 7, we show the distribution of disc mass versus local number density for the non-dispersed discs in all simulation runs. The left-hand panels show the discs at the end of the star formation process. The right-hand panels show the discs at the end of the simulations, |$2\, \textrm {Myr}$| after star formation has finished. The 2D distributions shown are the contours of the probability density, estimated using Gaussian kernel density on the logarithm of disc mass and stellar density. The 1D distributions are estimated using histograms of the logarithm of the respective quantity.
The distribution of local stellar density and gas mass of three disc populations, at two moments in the simulation. The left-hand panels show the population at the end of star formation, the right-hand panels at the end of the simulation. The populations are: discs born between |$0\, \textrm {Myr}$| and |$2.5\, \textrm {Myr}$| (purple, 6814 and 1645 discs at the end of star formation and the simulation, respectively); discs born between |$2.5\, \textrm {Myr}$| and |$3\, \textrm {Myr}$| (blue, 4465 and 1410 discs); and discs born between |$3\, \textrm {Myr}$| and |$6\, \textrm {Myr}$| (green, 6777 and 2754 discs).
The distribution of local stellar density and gas mass of three disc populations, at two moments in the simulation. The left-hand panels show the population at the end of star formation, the right-hand panels at the end of the simulation. The populations are: discs with radii between |$0\, \textrm {au}$| and |$20\, \textrm {au}$| (purple, 8216 and 1437 discs at the end of star formation and the simulation, respectively); discs with radii between |$20\, \textrm {au}$| and |$50\, \textrm {au}$| (blue, 6366 and 1260 discs); and discs with radii between |$50\, \textrm {au}$| and |$2000\, \textrm {au}$| (green, 3474 and 3112 discs).
Furthermore, we separate the discs in three populations, based moment of their birth (relative to the start of the simulation) in Fig. 6 and their disc radius in Fig. 7. These bins were chosen to have comparable numbers of discs in all bins at both times.
We calculate the local stellar density using the method by Casertano & Hut (1985) with the five nearest neighbours. While star formation is still ongoing, stars and discs form in regions spanning the whole range of stellar density. In particular, given our implementation of the star formation process with sink particles, many stars tend to form in regions of high stellar densities. Once the star formation process ends and the gas is expelled, the clusters expand. This results in a decrease in the overall density. The consequence is that discs are less often harassed dynamically, but also that photoevaporation becomes less effective. By the end of the simulations, almost no discs are present in regions of local density ≳ 105 stars pc−3.
In the bottom panels of Fig. 6 we can see how the local stellar density evolves. At the end of star formation, the stars born late (3–6 Myr) have not yet had time to evolve far from the gas distribution they formed in and form a rather symmetric, almost lognormal density distribution, similar to that of turbulent gas (Krumholz 2014). The older populations form a less symmetric distribution, with a heavier tail towards low densities. At the end of the simulation, each population has virialized and the density distributions of each of the populations are similar.
A similar trend is visible in the local stellar density of the disc radius bins. At the end of star formation, large discs are found in low-density regions, and small discs in high-density regions. 2 Myr later, the distributions of the three bins are similar.
At the end of star formation, the disc mass distributions of the three age bins peak at very similar values, but are subtly different at extreme values. The greatest disc masses belong to young discs, and smaller disc masses belong mostly to older discs. However, at the end of the simulation the distributions become more similar. The number of discs in each population differs considerably: of the old population about a quarter of discs survive, while of the young population about 40 per cent survives.
There is a clear trend for more massive discs to have larger radii. This is not simply a result of the initial conditions and the viscous expansion of the discs.
At the start of the simulations, the largest disc had a radius of 41 au. All discs that exceed this value at a later time have viscously grown to that size. External photoevaporation (and truncation) causes the discs to shrink, because we remove the evaporated gas from the outside. This process opposing viscous growth. At the end of our simulations the largest discs tend to be those that lose little mass and, therefore, are also among the most massive.
In Fig. 8, we present the median disc mass as a function of time. The solid lines show the times while star formation is still ongoing, and the dotted lines once the last star has formed. The mean disc mass varies during star formation because new discs are still forming and directly exposed to external photoevaporation. Once star formation stops, discs also stop forming, and, as a result, the median disc mass generally drops. In run #4, the median disc mass remains relatively constant, and in run #5, it even increases. This is possible if low-mass discs are preferentially dispersed while high-mass discs lose little mass.
Median disc mass in time for each simulation run. The solid lines correspond to ongoing star formation, and the dotted lines after it has ended. The shaded areas span the range between the 16th percentile and the 84th percentile. For clarity, these ranges are shown only for 2 runs, but the other ones span similar magnitudes.
In Fig. 9 we present the binned mean gas and solid masses of the discs as a function of the projected local stellar number density at the end of the simulations. The local stellar density is calculated in the same way as for Fig. 6, but with the distances between stars projected to two dimensions. The binned mean is calculated using a rolling bin spanning 100 stars. Both gas mass and solid mass are relatively constant across densities.
Binned mean gas mass (top panel, MJup) and solid mass (bottom panel, M⊕) at the end of the simulations versus projected local stellar density. The binned mean is calculated using a rolling bin spanning 100 stars. The local stellar number density is calculated as specified in Section 3.2.
In Fig. 10 we show the cumulative distribution of disc gas and solid masses at the end of star formation and the end of the simulations. Note that we normalize to all stars that have hosted discs. If the cumulative distribution has a maximum of 0.6, that means 60 per cent of stars retain discs. The lines show the cumulative distribution of discs from all our simulations. The shaded regions are between the extremes of the distributions of the individual simulations. By the end of star formation, almost all remaining discs have masses in solids ≳ 20M⊕, up to ∼600M⊕. After |$2\, \textrm {Myr}$| of evolution, the distribution spans almost the same range. All of the final discs have mass insolids in excess of 10 M⊕, which is a lower limit mass for the formation of rocky planets and the cores of gas giants (Ansdell et al. 2016). The masses obtained in the present simulations are higher than those in Concha-Ramírez et al. (2021), suggesting that the extended period of star formation may mediate the survival of massive discs. We expand on this discussion in Section 4.
Cumulative distribution of disc solid and gas masses at the end of star 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 the standard deviation.
4 DISCUSSION
4.1 Comparison with previous work
In previous work (Concha-Ramírez et al. 2019b, 2021), we performed simulations of star clusters in which stars with masses M* ≤ 1.9 M⊙ have circumstellar discs. These discs were subject to viscous evolution, external photoevaporation, and dynamical truncation. We assumed discs to be composed of 1 per cent dust (by mass). In those previous simulations, we found that discs become depleted of mass sufficiently quickly that |$\sim 60{{\ \rm per\ cent}}$| to 90 per cent are dispersed within |$2\, \textrm {Myr}$|. In these simulations, all stars were born instantaneously in a virialized Plummer distribution.
Here, we improve these models in two main ways: First, we start with a giant molecular cloud that forms stars through hydrodynamical collapse. This improved treatment of cluster formation leads to non-spherical and non-equilibrium initial stellar distributions. It also leads to stars being born over a time frame of 2 to 5 Myr. A second improvement is to the circumstellar disc model. We implement a model for external photoevaporation to remove the discs’ dust (Haworth et al. 2018a) and for internal photoevaporation to remove the discs’ gas (Owen et al. 2012; Picogna et al. 2019).
These improvements to the initial model by Concha-Ramírez et al. (2019b, 2021) resulted in interesting differences in the results. During the first few million years of evolution, the star formation process is ongoing and, as discs lose mass due to photoevaporation, the new discs that are constantly being formed keep the median mass of the overall population relatively constant (Fig. 8). As soon as this constant replenishing of discs stops, in most runs the median mass of the discs quickly drops (we discuss the runs with a different trend below). This is consistent with the results from Concha-Ramírez et al. (2019b, 2021). However, the mass in solids in our simulations (see the bottom panel of Fig. 9) is much higher than in Concha-Ramírez et al. (2021) in regions of comparable stellar projected density.
This discrepancy is partly due to a difference in IMF; our previous work used a lower limit of 0.01 M⊙ (which includes brown dwarfs), while here we use a lower limit of 0.08 M⊙ (the hydrogen burning limit). This explains why the average final solid masses in this work are even higher than the average initial dust masses in our previous work. However, the decrease in dust photoevaporation also plays a role, essentially locking in the solid content after ∼1 Myr. In Concha-Ramírez et al. (2021), the dust-to-gas ratio was assumed to be constant throughout the simulation.
4.2 Star formation recipe
The spatial distribution of star-forming regions also impacts disc-mass distributions. In previous work, we performed simulations starting from a spherically symmetric star cluster. Here we model the collapse of a molecular cloud to improve on these initial conditions. The resulting clusters, however, are, on average, smoother than observed star-forming regions (see Fig. 4). Stellar feedback, in particular stellar winds and jets from protostars, can substantially affect the morphology of star-forming regions (e.g. Mac Low & Klessen 2004; Hansen et al. 2012; Offner & Arce 2015).
Moreover, the SFE in our model is overestimated. While we set up an SFE of 30 per cent for the star formation process, the integrated SFE of molecular clouds may be orders of magnitude smaller (e.g. Chevance et al. 2020), and efficiencies of the order used in this work are found within much smaller clumps (∼0.1 pc radius, see e.g. Matzner & McKee 2000). The star formation efficiency per free fall time is also overestimated. The free fall time of our initial cloud is ∼1 Myr, and we form 30 per cent of the initial gas mass in stars over a few Myr. This implies an efficiency per free fall time of ∼10 per cent. However, measurements of this quantity through various traces indicate a value closer to 1 per cent (Krumholz et al. 2019).
For this reason, we ran two additional series of simulations where we varied the recipe for star formation, but used the same initial conditions and realizations of the IMF.
In the first series, we lowered the absolute star formation efficiency from 0.3 to 0.1. These runs were identical to the runs presented above (up to propagating numerical errors) until they reached the prescribed star formation efficiency, between 1 and 3 Myr before the runs with higher efficiency. We refer to these runs as the low-efficiency runs.
We present the distribution of the discs’ local stellar density and gas mass in Fig. 11. We compare three populations at the same moment in time, namely the end of the low-efficiency runs, 2 Myr after the end of the star formation process. The three populations are the discs in the low-efficiency runs, the discs in the high-efficiency runs born before the moment star formation ended in the low-efficiency runs (effectively the equivalent of the low-efficiency runs’ populations), and the all discs in the high-efficiency runs (born until that moment; star formation has not necessarily ended in these runs).
The distribution of local stellar density and gas mass of three disc populations. The populations are: all discs in the low-efficiency runs (purple, 4064 discs); the discs in the high-efficiency runs born early (blue, 3598 discs); and all discs in the high-efficiency runs (green, 15130 discs). The data from the low-efficiency runs are from the end of the runs. The data from the high-efficiency runs are from the moment that the low-efficiency run with the same initial conditions ends.
Comparing the full populations of both sets of runs, we can see that the disc mass distributions peak around the same mass (∼10 MJup), but that of the high-efficiency runs have more discs at higher disc masses, while the low-efficiency runs have more discs slightly below the peak. The low mass tails are indistinguishable. Between the population of the low-efficiency runs and the first discs of the high-efficiency runs, the disc–mass distribution of the low-efficiency populations peaks at slightly higher masses than the early high-efficiency populations. Two effects compete here: The first is the late formation of discs in the high-efficiency run, which leads to the presence of more massive discs (corresponding to those formed recently). The second is the formation of a larger number of massive stars in the high-efficiency runs. This leads to more efficient external photoevaporation of the early population of the high-efficiency runs. This is visible in the slight shift of the peak of the disc mass distribution.
In the second series of runs, we also lowered the rate of star formation, by increasing the parameter td in equation (1) to 10 Myr. We kept the star formation efficiency at 0.1. These runs differ from the main runs presented from the moment stars form. We refer to these runs as the low-rate runs. In the runs with a low star-formation rate, the epoch of star formation was a few Myr longer than in the runs with high star formation. Together with the decreased absolute SFE, this puts the SFE per freefall time close to the observed value.
We present the distribution of the discs’ local stellar density and gas mass in Fig. 12. We compare three populations: the discs in the low-efficiency runs at high and low rates, at the moment the high-rate runs end; and those in the low-efficiency, low-rate runs, at the end of those simulations.
The distribution of local stellar density and gas mass of three disc populations. The populations are: discs in the low-efficiency, high-rate runs at the end of those simulations, which we call t1 (purple, 4064 discs); discs in the low-efficiency, low-rate runs at t1 (blue, 4458 discs); and discs in the low-efficiency, low-rate runs at the end of those simulations, which we call t2 (green, 1579 discs).
The disc mass distribution of the low-rate runs compared to the high-rate runs at the same moment is shifted to greater masses. At the end of the low-rate runs, the distribution has shifted towards lower masses. Again, it appears that a population with more young discs also has more massive discs.
4.3 Comparison with observations
Ansdell et al. (2016) propose that a rapid depletion of gas mass in discs can lead to the observed characteristics of exoplanet populations. Traditional theories of planet formation indicate that ∼10M⊕ cores should be able to accrete gaseous envelopes at a rate of ∼1MJup within |$\sim 0.1\, \textrm {Myr}$|. Observational surveys, on the other hand, indicate that ‘super-Earths’, or intermediate-mass rocky planets, are about an order of magnitude more common than gas giants (e.g. Petigura, Howard & Marcy 2013). Ansdell et al. (2016) argue that, if typical ∼10M⊕ cores form in discs that are already depleted of gas, this contributes to the observed diversity in the number of planetary types. Results by Lopez & Rice (2018) indicate that late (≳ |$10\, \textrm {Myr}$|) rocky planet formation could explain the newly observed population of highly irradiated rocky planets of ∼1.5R⊕, supporting the idea that rocky planets form in discs which have already lost most of their gas. The evolution of dust and gas disc masses in our simulations suggests a similar context for planet formation. Sellek, Booth & Clarke (2020) demonstrate that photoevaporation of dust leads to a decrease in pebble flux in the inner disc after the pebble flux has already peaked, but before 1 Myr. This implies that external photoevaporation of dust need not inhibit the formation of cores.
We do note that our simulations are of massive, high-density star forming regions with massive stars. Locally, the various Orion star forming regions are analogues. However, other regions such as the Lupus and Taurus clouds are of lower mass and density, and lack massive stars. The FUV radiation field in such regions is typically small, and may even be dominated by stars outside the region. Cleeves et al. (2016) estimate that one member of the Lupus clouds may be receiving a radiation field of ∼3 G0 due to nearby OB stars. Haworth et al. (2017) then showed that this radiation field is driving a photoevaporating flow. The photoevaporation rates were computed using a method predating the more carefully modelled FRIED grid. Because the FRIED only provides photoevaporation rates for radiation fields down to 10 G0, we are unable to self-consistently model low-mass star forming regions.
Circumstellar discs are observed in regions where external radiation should have evaporated them already. This ‘proplyd lifetime problem’ is notoriously present in the ONC, where discs observed within |$0.3\, \textrm {pc}$| of the massive star θ1C Ori are exposed to strong radiation fields (∼104 G0, e.g. O’dell & Wen 1994). The fact that there are still discs observed in the vicinity of the star suggest that either the discs were initially massive (≳ 1 M⊙), or that θ1C Ori is younger than |$0.1\, \textrm {Myr}$|. Given that such massive discs are gravitationally unstable and that the stellar age distribution in the ONC is |$\sim 1\, \textrm {Myr}$| (e.g. Da Rio et al. 2010, 2012), there must be other factors at play for these discs to have survived. Störzer & Hollenbach (1999) argue that these discs are observed at a special moment in time: while coincidentally passing close to the centre of the ONC, but they have spent the majority of their lives in regions of lower background radiation fields. Scally & Clarke (2001) simulate the inner area of the ONC and propose that the type of orbits necessary for the hypothesis of Störzer & Hollenbach (1999) to be feasible are not dynamically possible. However, a short-lived trajectory through high stellar density regions might not be the only way to explain the presence of massive discs in the centre of the ONC. Winter et al. (2019) propose a combination of mechanisms, among which an extended period of star formation, that enable the presence of massive discs in the vicinity of massive stars, even at relatively late times. Our results also show that the high mass component of a disc population typically consists of young discs.
Recently, several studies (Kruijssen, Longmore & Chevance 2020; Winter et al. 2020b; Chevance, Kruijssen & Longmore 2021; Longmore, Chevance & Kruijssen 2021) reported correlations between the architecture of exoplanetary systems and the phase space densities in the vicinity of their host stars. These results imply that external influences can shape exoplanetary systems. The authors propose two origins of these differences in phase space density: traces of the star’s phase space density within its birth region, or large-scale galactic perturbations. These imply that exoplanetary systems are shaped by processes in the birth cluster and in the galaxy, respectively. Subsequent work by Kruijssen et al. (2021) showed that the phase space overdensities coincide with known structures driven by large-scale galactic perturbations. This favours the galactic perturbation scenario. The existence of influences on planet formation due to cluster processes are not ruled out, but the galactic population of field stars is sufficiently well-mixed to erase signatures of the structure of the star-forming region. Also, Kruijssen et al. (2012) proposes that bound clusters arise from star formation in regions of high ISM and stellar density, which implies that the unbound population (which forms the field population) consists of stars formed in low-density regions. The impact of a high-density birth environment on planet formation can perhaps be uncovered in the contrast between the exoplanetary systems of cluster stars and field stars.
5 SUMMARY AND CONCLUSIONS
We perform simulations of molecular cloud collapse and star formation. Our calculations begin with the collapse of a 104 M⊙ molecular cloud with a 3 pc radius. The star formation process ends when the star formation efficiency reaches 30 per cent. The excess gas is then instantaneously removed. From this point, we continue the simulations using a combination of N-body dynamics, dynamical truncation, viscous disc evolution, and stellar evolution. After formation stars with a mass M* ≤ 1.9M⊙ receive circumstellar discs. Stars more massive than M* > 1.9M⊙ are emitting UV radiation, which affects the discs in their vicinity. We also take the internal photoevaporation of the host star into account, as well as the dynamical truncations and viscous evolution. We also implemented a model for photoevaporation of dust. We run the simulations for |$2\, \textrm {Myr}$| after the last star has formed. We conclude that:
A prolonged period of star formation, lasting from 2 to |$5\, \textrm {Myr}$|, allows for relatively massive (Mgas ∼ 100 MJup, Mdust ∼ 100 M⊕) discs to survive the high intensity UV radiation of nearby cluster members, and therewith provide a solution to the proplyd lifetime problem. An extended periods of star formation then mediates the survival of massive discs in regions with a strong radiation field.
The discs that make it to the end of our simulations are preferentially the ones that were born later (after |$\sim 3\, \textrm {Myr}$|) in the star formation process.
While photoevaporation removes gas from the discs, dust quickly becomes resistant to photoevaporation. This might allow discs to keep enough mass to form rocky planets, even when depleted of gas.
Photoevaporation due to the radiation of nearby massive stars is an essential process that drives the dissolution of circumstellar discs. The structure of the star cluster is then important in determining the effectiveness of this process. Intracluster gas can shield discs from this process, in which case stellar dynamical processes become relatively more important. Another important aspect of cluster formation is the local star-formation history. Massive stars born late may be unable to evaporate circumstellar discs effectively, because the cluster, by that time, has expanded, and the discs may already have turned into planetary systems. On the other hand, a late-formed disc is protected by the larger distance to massive stars of the older cluster. Such lower density may be caused by the expulsion of the primordial gas from the parent cluster.
These relatively fundamental processes in the star-forming environment give rise to a wide variety of disc masses and sizes that depend on the dynamical history of the parent star. Two identical stars with identical discs born in the same cluster but on a different orbit or born at a different time can then develop completely different disc parameters. The resulting planetary system will subsequently also be very different. In our simulations we observe ranges in disc masses and sizes extending over more than three orders of magnitude, just from the time and orbit in which the disc-hosting star was born.
ACKNOWLEDGEMENTS
It is a pleasure to thank the referee for detailed comments. This project was supported by NOVA, and we used the National Dutch Supercomputer Cartesius, and Little Green Machine-II.
Software: The present works makes use of the following software: amuse (Pelupessy et al. 2013; Portegies Zwart et al. 2013), fi (Pelupessy et al. 2004), ph4, bridge (Fujii et al. 2007; Portegies Zwart et al. 2020), seba (Portegies Zwart & Verbunt 1996; Toonen et al. 2020), vader (Krumholz & Forbes 2015), numpy (Van Der Walt, Colbert & Varoquaux 2011), scipy (Virtanen et al. 2020), matplotlib (Hunter 2007), and makecite (Price-Whelan, Mechev & jumeroag 2018).
Energy consumption: The series of simulations presented in this work took ∼20 d to run on 30 cores each. Although we present the results of six runs, triple that number was eventually run (including code development, testruns, and referee requests for additional simulations). The total amount of computer time used then tops 259.200 h. Considering 12 Wh for a CPU, this results in ∼3000 kWh. Using the conversion factor 0.283 kWh kg−1 (Portegies Zwart 2020), would results in ∼10 000 kg of CO2. We used two computers for the calculations, the Dutch National supercomputer Cartesius, and our local 192-core workstation. Both machines are powered with renewable resources; the former using Norwegian hydroelectric power (using certificates) and the latter by Dutch wind-mill energy.
DATA AVAILABILITY
The code underlying the simulations in this article is available at https://doi.org/10.5281/zenodo.6579534.
