Effects of stellar density on the photoevaporation of circumstellar discs

ABSTRACT

Circumstellar discs are the precursors of planetary systems and develop shortly after their host star has formed. In their early stages, these discs are immersed in an environment rich in gas and neighbouring stars, which can be hostile for their survival. There are several environmental processes that affect the evolution of circumstellar discs, and external photoevaporation is arguably one of the most important ones. Theoretical and observational evidence point to circumstellar discs losing mass quickly when in the vicinity of massive, bright stars. In this work, we simulate circumstellar discs in clustered environments in a range of stellar densities, where the photoevaporation mass-loss process is resolved simultaneously with the stellar dynamics, stellar evolution, and the viscous evolution of the discs. Our results indicate that external photoevaporation is efficient in depleting disc masses and that the degree of its effect is related to stellar density. We find that a local stellar density lower than 100 stars pc⁻² is necessary for discs massive enough to form planets to survive for 2.0 Myr. There is an order of magnitude difference in the disc masses in regions of projected density 100 versus 10⁴ stars pc⁻². We compare our results to observations of the Lupus clouds, the Orion Nebula Cluster, the Orion Molecular Cloud-2, Taurus, and NGC 2024, and find that the trends observed between region density and disc masses are similar to those in our simulations.

1 INTRODUCTION

Circumstellar discs are the reservoirs of gas and dust that surround young stars and have the potential to become planetary systems. Their evolution will determine the time and material available to form planets. Studying the evolution of circumstellar discs can then help us understand planet formation and the diversity of observed planetary systems.

These circumstellar discs develop almost immediately after star formation, as a direct consequence of the collapse of a molecular cloud and angular momentum conservation (Williams & Cieza 2011). Their surroundings are rich in gas and neighbouring stars, which can be hostile to the discs and affect their evolution in different ways. In environments with high stellar densities, dynamical encounters with nearby stars can truncate the discs (e.g. Pfalzner, Olczak & Eckart 2006; Vincke, Breslau & Pfalzner 2015; Portegies Zwart 2016; Bhandare & Pfalzner 2019). Face-on accretion of gas on to the circumstellar discs can cause them to shrink and increase their surface densities (e.g. Wijnen et al. 2016, 2017). Feedback from processes related to stellar evolution, such as stellar winds and supernovae explosions, can also truncate, tilt, or completely destroy the discs (Pelupessy & Portegies Zwart 2012; Close & Pittard 2017; Portegies Zwart et al. 2018). The presence of bright, massive stars in the vicinity of circumstellar discs can heat their surface enough to evaporate mass from them. This process, known as external photoevaporation, is arguably one of the most important environmental mechanisms in depleting mass from young circumstellar discs, and its effects seem to greatly outperform that of other means for disc truncation (e.g. Hollenbach, Yorke & Johnstone 2000; Adams et al. 2004; Guarcello et al. 2016; Facchini, Clarke & Bisbas 2016; Winter et al. 2018; Haworth & Clarke 2019; Winter et al. 2020b; Haworth & Owen 2020).

The effects of external photoevaporation have been identified in observational surveys of young stellar objects in star-forming regions. Proplyds – cometary tail-like structures formed by ionized, evaporating discs – have been observed in particular in dense regions of the Orion nebula (O’dell & Wen 1994; O’dell 1998; Vicente & Alves 2005; Eisner & Carpenter 2006; Mann et al. 2014; Kim et al. 2016). Surveys of protoplanetary discs in star-forming regions show that discs closer to bright stars are less massive than their counterparts in sparser regions (Mann & Williams 2009; Fang et al. 2012; Mann et al. 2014; Ansdell et al. 2017; van Terwisga et al. 2020), suggesting that discs in the vicinity of these stars are strongly affected by their environment. Disc fractions (the number of young stellar objects around which dust is detected, over the total number of objects) and disc mass distributions in younger and less dense star-forming regions, such as Lupus and Taurus, are statistically indistinguishable from each other in terms of disc mass distributions. The average disc mass in these regions is higher than in the Orion Nebula Cluster (ONC, Ansdell et al. 2016; Eisner et al. 2008, 2018; van Terwisga, Hacar & van Dishoeck 2019), which is a much denser environment.

The extent of the effects of radiation in depleting disc mass depends on the proximity to bright stars. As such, the effects of external photoevaporation are important in clustered environments, and different theoretical models have been developed to study the process in that context. Scally & Clarke (2001) model the dynamics of a cluster with 4000 stars, with discs of radii 100 au which remain constant throughout the simulation. The mass loss due to photoevaporation is calculated in post processing, by keeping track of the radiation to which each star is exposed during the simulation. Their results show that external photoevaporation is important in depleting disc mass in regions similar to the centre of the ONC. A different approach is taken by Adams et al. (2006) and Fatuzzo & Adams (2008), who model the dynamics of star clusters of different sizes and derive the background far-ultraviolet (FUV) radiation for each simulation. They then estimate the photoevaporation mass-loss rates of the discs depending on the average background radiation that hey have been exposed to. Adams et al. (2006) model star clusters with 100, 300, and 1000 stars and find that external photoevaporation is only important for disc radii larger than 30 au, due to the low average background UV exposure. Furthermore, models of single, externally illuminated discs show that the supersonic flows caused by FUV photons heating the disc surface can explain the observed proplyd shapes (e.g. Richling & Yorke 1997; Johnstone, Hollenbach & Bally 1998; Störzer & Hollenbach 1999; Adams et al. 2004).

In more recent work, Winter et al. (2018, 2020a) use a statistical approach to model background FUV radiation fields in regions of different stellar number density. Winter et al. (2020a) find that 90 per cent of circumstellar discs are destroyed by external photoevaporation within 1.0 Myr in a region comparable to the Central Molecular Zone of the Milky Way, and that the effects of photoevaporation are particularly destructive for discs around low-mass stars (M_* < 0.3M_⊙). For regions similar to the solar neighbourhood (surface density Σ₀ = 12 M_⊙ pc⁻²) they find a mean dispersal time-scale of |$\sim {3.0}\, {\rm Myr}$|⁠. Similar results are obtained by Nicholson et al. (2019), who calculate mass-loss rates in N-body simulations using the same prescriptions as Scally & Clarke (2001) and find that, in regions with high degree of substructure (density ∼100 M_⊙ pc⁻³), 50 per cent of the discs with initial radii ≥100 au are destroyed by external photoevaporation within 1 Myr. In regions of lower densities (∼10 M_⊙ pc⁻³), half of the discs are destroyed within 2 Myr.

In Concha-Ramírez et al. (2019, hereafter Paper I), we presented a new numerical implementation which allows for self-contained simulations of external photoevaporation in clustered environments. Photoevaporation process is solved simultaneously with the stellar dynamics (including disc encounters and truncations), stellar evolution, and viscous spreading of the circumstellar discs. This causes disc masses, sizes, and column densities to vary in time, and the mass-loss rate of the discs is calculated accordingly. The results of the simulations in Paper I show that external photoevaporation is efficient in destroying circumstellar discs on a relatively short time-scale. For regions of stellar densities ∼100M_⊙pc⁻³, around |$80{{\ \rm per\ cent}}$| of discs have evaporated within 2.0 Myr. Between 25 per cent and 60 per cent of the discs, depending on region density, are destroyed within the first 0.1 Myr. We argue that the rapid decrease in disc mass is dominated by external photoevaporation, rather than dynamical truncations, and that the former mechanism constrains the time available for planet formation.

Observational and theoretical evidence suggest that the local stellar density is a key factor in the survival of circumstellar discs and in their eventual observed mass distributions. Understanding disc mass and size distributions in young star clusters is therefore paramount for understanding planet formation and evolution. Here, we use the model developed in Paper I to determine for which range of stellar densities the effects of photoevaporation are most efficient. We perform simulations of circumstellar discs embedded in star clusters and explore a parameter space of stellar densities spanning five orders of magnitude. The clusters are evolved for 2.0 Myr and we investigate the final mass distributions of the disc population. We compare our simulation results to observed dust masses of young stellar objects in the Lupus clouds (Ansdell et al. 2016, 2018), the ONC (Mann et al. 2014; Eisner et al. 2018), the Orion Molecular Cloud-2 (van Terwisga, Hacar & van Dishoeck 2019), the Taurus region (Andrews et al. 2013), and NGC 2024 (Getman, Feigelson & Kuhn 2014; van Terwisga et al. 2020).

2 MODEL

We use the Astrophysical Multipurpose Software Environment (amuse¹, Portegies Zwart et al. 2009, 2013), to bring together codes for viscous disc evolution, stellar dynamics, and stellar evolution, along with an implementation of external photoevaporation. The setup and models used for the simulations in this paper are the same as in Paper I. In this work, we perform simulations spanning a larger range of stellar densities.

Here, we present a summary of the implementation used for the simulations. For a detailed explanation of the model the reader should refer to Paper I. All the code developed for the simulations, data analyses, and figures of this paper is available online.²

2.1 Stars and circumstellar discs

We separate the stars in the simulations into two populations: stars with masses M_* ≤ 1.9 M_⊙, and stars with masses M_* > 1.9 M_⊙. The reason for this mass limit is related to the photoevaporation mass-loss calculation and further explained in Section 2.2. This mass separation is for photoevaporation purposes only and does not influence the dynamical evolution of the stars. All stars with masses M_* ≤ 1.9 M_⊙ are surrounded by a circumstellar disc, while stars with higher masses have no discs and are considered only as generating ionizing radiation. Massive stars are subject to stellar evolution, implemented using the code seba (Portegies Zwart & Verbunt 1996; Toonen, Nelemans & Zwart 2012) through its amuse interface. Stars with discs do not undergo stellar evolution in the simulations. The dynamical evolution of the clusters is implemented using the fourth-order N-body code ph4, incorporated in amuse.

Circumstellar discs are implemented using the Viscous Accretion disc Evolution Resource (vader) developed by Krumholz & Forbes (2015). vader models mass and angular momentum transport on a thin, axisymmetric disc. This allows us to take into consideration the viscous spreading of the discs. Each vader disc in our simulations is composed of a grid of 100 logarithmically spaced cells between 0.05 and 2000 au. The discs have a turbulence parameter of α = 5 × 10⁻³.

2.2 External photoevaporation

External photoevaporation is dominated by FUV photons (Armitage 2000; Adams et al. 2004; Gorti & Hollenbach 2009). To model the FUV radiation from the massive stars we pre-compute a relation between stellar mass and FUV luminosity using the uvblue spectral library (Rodriguez-Merino et al. 2005). The obtained FUV luminosity fit is shown in fig. 2 of Paper I. During the simulations, we use this fit to obtain the FUV luminosity of each massive star at every time-step.

Mass loss due to external photoevaporation is calculated for each disc using the fried grid (Haworth et al. 2018b). The fried grid provides a set of pre-calculated, external photoevaporation mass-loss rates for discs immersed in radiation fields of varying intensity, from 10 to 10⁴ G₀. The grid spans discs of mass ∼10⁻⁴ to 10² M_Jup, radius from 1 to 400 au, and host star mass from 0.05 to 1.9 M_⊙. To stay within the limits of the grid, we give all stars with masses M_* ≤ 1.9 M_⊙,  a circumstellar disc, and all stars with masses M_* > 1.9 M_⊙ are considered as only generating radiation.

We calculate the mass loss of every disc as follows. For each disc we begin by calculating its distance to every star of mass M_* > 1.9 M_⊙ and determining the total radiation that the disc receives from those stars. We do not consider extinction in this calculation. We then use this total radiation and the disc parameters to interpolate a mass-loss rate |$\dot{\mathrm{M}}$| from the fried grid. This |$\dot{\mathrm{M}}$| is then used to calculate the total mass lost by the disc in the current time-step. Assuming a constant mass-loss rate over the time-step, the mass is removed from the outer regions of the disc: we advance over the disc cells starting from the outside removing mass from each, until the corresponding amount of mass has been removed. Through this process, mass loss due to photoevaporation results in a decrease of disc mass and disc radius.

We consider a disc dispersed when its mass drops below 0.03 M_⊕, based on the non-detection mass limits from (Ansdell et al. 2016), or when its mean column density is lower than 1 g cm⁻² (Pascucci et al. 2016, see also Section 2.4).

2.3 Initial conditions

2.3.1 Star clusters

We simulate clusters with 10³ stars and initial virial radii of 0.1, 0.3, 0.5, 1.0, 2.5, and 5.0 pc. Stars are initially distributed in a Plummer (1911) sphere. Stellar masses are drawn from a random Kroupa (2001) mass distribution with upper limit 100 M_⊙. All models start in virial equilibrium (virial ratio Q = 0.5). No primordial mass segregation, binaries, or higher multiplicity systems are considered.

In Table 1, we present the models used for this work. The mean number of stars with discs in each simulation is 974.7 ± 1.7. The mean mass of the stars with discs is |$0.23\substack{+1.66 \\-0.22}\,\mathrm{M}_{\odot }$|⁠. The mean number of stars generating UV radiation is 25.3 ± 1.7. The third column of Table 1 shows the mass ranges spanned by these stars.

We evolve each cluster for 2.0 Myr. We run each model six times, with a different random seed for the mass function and the initial stellar positions and velocities.

2.3.2 Circumstellar discs

2.4 Model caveats

Our model is designed as a controlled experiment to investigate the physical processes going on inside star-forming regions, in particular with regards to external photoevaporation. There are quite a number of assumptions in our simulations which we justify based on previous theoretical work and observations. Below we discuss some of these processes and parameters and the implications they might have.

Star-forming regions are not only rich in stars but also in gas, which can linger for several million years (Portegies Zwart, McMillan & Gieles 2010). Intracluster gas could influence our results in two main ways: first, the presence of gas and its subsequent expulsion in time affect the virial equilibrium and thus the dynamics of the star clusters. Second, gas can absorb some of the FUV radiation coming from bright stars, effectively protecting the discs from external photoevaporation and allowing them to live for longer (Winter et al. 2019; Ali & Harries 2019; van Terwisga et al. 2020), therefore giving more time for the planet formation process to occur. The presence of gas can also explain the ‘proplyd lifetime problem’ observed in the ONC, in which discs not massive enough to survive in the environment of θ¹C Ori are still observed in the region (Winter et al. 2019).

The Far-ultraviolet Radiation Induced Evaporation of Discs (fried) grid was constructed using a 1D disc model, but later simulations by Haworth & Clarke (2019) show that mass-loss rates can increase up to a factor of 3.7 when considering 2D discs. It is likely then that the mass losses used in this work are only a lower limit for the effects of external photoevaporation.

Internal photoevaporation, the process in which X-ray and UV photons coming from the host star itself lead to mass loss, is not considered in these simulations. Internal photoevaporation can drive mass loss in the inner regions of the discs (⁠|$\sim 1 - {10}$| and 30 au, Gorti, Dullemond & Hollenbach 2009; Gorti & Hollenbach 2009) and even in outer regions under certain conditions (Font et al. 2004; Owen et al. 2010). However, external photoevaporation is arguably the dominant process in regions |$\gt {10}\, {\rm au}$| (Hollenbach, Yorke & Johnstone 2000; Fatuzzo & Adams 2008). Our approximation of external photoevaporation removing mass from the outer regions of the disc only is also idealized, since while internal photoevaporation seems to clear the disc at specific radii, FUV photons coming from external sources can heat and evaporate mass from the whole disc surface (Adams et al. 2004).

Our conditions for disc dispersal might be overestimating the number of destroyed discs, particularly for relatively high-mass stars. In some cases, a disc of density 1 g cm⁻² or lower can still be detectable with modern instruments. As the mean mass of stars with discs in our simulations is 0.23 M_⊙ (corresponding to an initial disc mass of 24 M_Jup), and as most of the massive discs in our simulations survive until the end (see Section 3.1), we consider this possible underestimation of surviving discs to be within the uncertainty of our simulations and to not affect our results.

Our prescription of the EUV mass-loss rate (equation 2) is taken from Johnstone et al. (1998) and corresponds to a thick photodissociation region (PDR). In practice, in EUV-dominated fluxes the PDR is expected to be thin. This can change the value of the EUV mass-loss rate |$\dot{M}_{\mathrm{ EUV}}$| calculated in each time-step. However, since our focus is on FUV photoevaporation, and since discs are only in the EUV regimes for a short time in our simulation, this hardly affects our mass-loss rate calculations.

Regarding disc masses, it is generally accepted that a 100:1 gas-to-dust mass ratio defines the composition of circumstellar discs. However, several authors have pointed out that this value might change across discs and in time (Williams & Best 2014; Manara et al. 2020). This can lead to observed disc dust masses being greatly underestimated (Manara, Morbidelli & Guillot 2018). New models of externally irradiated, evaporating discs by Haworth et al. (2018a) show that considering grain growth can lead to less dust being lost through external photoevaporation, and thus to the dust:gas ratio increasing in time. A more careful implementation of the separate dust and gas components in a disc can help to overcome this problem.

The distribution of stars in a Plummer sphere is an idealized geometry. Star-forming regions have complex configurations and can present fractal structures, filaments, and other regions of increased surface density (e.g. Scalo 1990; Elmegreen & Falgarone 1996; Elmegreen et al. 2000; Bate 2010; Hacar et al. 2013; Chevance et al. 2020; Krause et al. 2020). The simulations carried out for this work represent only local densities, but for improved analyses of disc survival in star-forming regions it is important to consider different spatial distributions.

3 RESULTS

To illustrate the evolution of our model in time for individual discs, in Fig. 1 we show the evolution of several stars and their corresponding circumstellar discs. These particular tracks are taken from one of the realizations of model R1.0 (see Table 1). We show seven stars with discs as they move through the cluster. Black crosses mark the position of each star at the beginning of the simulation, and the label next to each shows the mass of the star. The sizes of the coloured circles in the stellar tracks are proportional to the disc radii, and their colour indicates the total disc mass. Red crosses, where present, show the moment when the disc is dispersed. The black thin lines that follow a red cross indicate the continuation of the orbit of the star, which keeps moving through the cluster after its disc has been evaporated. The trajectories of some massive, radiating stars are shown in thin blue lines. The solid and dashed circles in the background show the core radius and half-mass radius of the cluster, respectively, at |$t = {0.0}\, {\rm Myr}$|⁠. A disc around a 0.37 M_⊙ star survives all through the simulation, however, the variations of its radius in time due to photoevaporation and viscous expansion can be seen. A 0.14 M_⊙ star initially near the centre keeps its disc until around halfway through the simulation. A very low-mass star, 0.06 M_⊙, loses its disc very quickly even if located in the periphery of the cluster. While our simulations are 3D, in this illustrative figure we show a 2D projection of the location of the stars.

3.1 Disc fractions and lifetimes

We define the disc fraction at time t as the number of discs at t over the initial number of discs in each cluster. In Fig. 2, we show disc fractions separated in terms of the mass of their host stars: low-mass stars (M_* ≤ 0.5 M_⊙) in the top panel and high-mass stars (0.5 M_⊙ < M_* ≤ 1.9 M_⊙) in the bottom panel. The disc fraction for high-mass stars stays constant through time for the R1.0, R2.5, and R5.0 models. These discs lose mass but not enough to be completely evaporated, except for a slight decrease near the end for the R1.0 model. In the R0.1, R0.3, and R0.5 models, however, starting around 1.0 Myr even massive discs get destroyed. Final disc fractions decrease with increasing stellar number density. The R0.1 and R0.3 models show a very similar evolution, meaning that the density of the R0.3 model is a higher limit for the effects of external photoevaporation in destroying discs in such simulations.

A large drop in the number of stars with discs before 0.2 Myr is observed in models R0.1 and R0.3. Similar behaviour was obtained in the simulations performed in Paper I. This drop is related to discs around very low-mass stars being dispersed rapidly once the simulation begins and external photoevaporation is ‘switched on’. This drop can be seen in all the curves, but it becomes less pronounced for lower densities.

In Table 2, we present the mean disc lifetimes and disc half-life for each model, averaged over six runs. Disc lifetimes are calculated as a mean of the times when a disc is completely dispersed in the simulations, following the dispersion criteria explained in Section 2.2. It is important to mention that this mean is calculated only considering the discs that get dispersed within the 2.0 Myr spanned by the simulation, and the discs that survive would likely increase these values. Considering the resulting disc fractions, however, the obtained disc lifetimes are a good approximation. The disc half-life corresponds to the moment when half of the discs in a simulation run have been dispersed. Both the disc lifetimes and the half-life increase with decreasing stellar density.

3.2 Disc masses

In Fig. 3, we show the evolution of the mean disc mass in time versus the local stellar number density. The local stellar number density is calculated for each star using the method described by Casertano & Hut (1985) with the five nearest neighbours. The binned mean disc mass is calculated using a sliding bin spanning 100 stars.

The thick dotted lines in Fig. 3 show the mean disc mass at |$t = {0.0}\, {\rm Myr}$|⁠, and the thick solid lines at t = 2.0 Myr. The shaded areas around these lines represent the standard error. The thin lines in between show the evolution of the curve in 0.2 Myr intervals. The expansion of the clusters in time is reflected by the t = 2.0 Myr curves spanning larger density ranges than the |$t = {0.0}\, {\rm Myr}$| curves. This effect is less pronounced in the R1.0, R2.5, and R5.0 models because they expand in a longer time-scale. The slope of the final mean disc mass distribution increases with decreasing stellar density. This is related to the core density in each region, which is also decreasing: the curves in the R5.0 model are several orders of magnitude lower, in terms of density, than the R0.1 model. The R0.1 model has a distribution of disc masses such that the most massive discs are found further away from the centre, with differences of about one order of magnitude between the discs located in the centre and in the outskirts of the cluster. In the R5.0 model, the mass difference between discs in different locations is much smaller, and the disc masses are of the same order of magnitude through all the density range.

In Fig. 4, we show the mean dust mass of the discs versus the projected local stellar density. We use a 1:100 dust:gas mass ratio to determine the dust mass of our discs. We calculate the density in the same way as in Fig. 3, but projecting the distances between stars to two dimensions. This allows us to compare disc dust masses in our simulations to observed disc populations. The blue diamonds show points of observed average disc dust mass versus local density of young stellar objects for the Lupus clouds (Ansdell et al. 2016, 2018), the ONC (Mann et al. 2014; Eisner et al. 2018), the Orion Molecular Cloud-2 (OMC-2, van Terwisga et al. 2019), Taurus (Andrews et al. 2013), and NGC 2024 East and West (Getman, Feigelson & Kuhn 2014; van Terwisga et al. 2020), as labelled.

Fig. 4 shows a break in disc masses around a local density of 100 stars pc⁻². The slope of the disc mass distribution changes around that point for all models, except R5.0. In models R2.5 and R1.0 we see the slope of the distributions increasing as we move towards higher densities. For the R0.1, R0.3, and R0.5 models, we see disc mass distributions stay relatively constant for densities lower than 100 stars pc⁻², and for higher densities we see a negative slope. The difference in masses for discs in regions of density between 100 and 10⁴ stars pc⁻² is about one order of magnitude. A similar effect can be seen in the observational points, except for NGC 2024 East. This behaviour suggests that 100 stars pc⁻² is a critical density for determining disc masses.

The average disc dust masses of the observations are calculated by fitting a lognormal distribution on the masses. Although disc masses span a large dynamic range, a lognormal distribution is a good description of these populations (Williams et al. 2019). The local density for each point is calculated using the five nearest neighbours method. Lupus data are an average for all the clouds, using the complete list of Class II sources in Ansdell et al. (2016, 2018). It is important to note that the Lupus III cloud dominates the signal for that particular region, because it has the largest population of Class II sources. For the OMC-2, the data come from van Terwisga et al. (2019), who use the source catalogue from Megeath et al. (2012) assuming completeness. ONC data come from Megeath et al. (2016), including completeness corrections.

In the OMC-2 and ONC regions, observations sample two different density regimes in the same cloud, relatively close together in space. Therefore, the conditions in our models most closely resemble the properties of the discs that were sampled by the observations, and we can interpret them as different density bins along a single model. It is immediately apparent that both the gradient of average disc mass with density as well as the average disc masses themselves resemble the models closely. Given the considerable uncertainties in extracting disc masses from millimetre-continuum observations (see, for instance the discussion in Eisner et al. 2018) the similarity in the gradients suggests that our models are successful at capturing the general behaviour of external photoevaporation.

In NGC 2024, Getman et al. (2014) find evidence for an age gradient of young stars, which van Terwisga et al. (2020) suggest as an explanation for the large difference in mean disc masses in NGC 2024 East and NGC 2024 West. In NGC 2024, the western part is the older and resembles the ONC in age and conditions, while the eastern disc population has a lower average age. We represent this difference in the figure by making NGC 2024 East in a different shade. Comparing these models to the observations, we see that the NGC 2024 West data lie closely to those of the ONC, while the data for the eastern subpopulation occupy a region of average disc mass and local stellar density space which is more consistent with a younger, compact population of stars.

Lupus and Taurus discs, on the other hand, sample a much more heterogeneous part of parameter space in terms of initial densities. Our models do not closely resemble the conditions under which the stars in these samples formed (see, for instance Roccatagliata et al. 2020). However, the result that when the average stellar number density is low enough (below ∼100 stars pc⁻²), the average disc masses are similar at similar ages does seem to apply to these star-forming regions, even though this is a part of parameter space we do not explore.

4 DISCUSSION

In the simulations performed for this work, external photoevaporation is effective in destroying the majority number of discs within 2.0 Myr of evolution. The initial stellar density of each region affects the fraction of surviving discs, as well as their final mass distributions. In all models, except for R5.0, half of the discs are destroyed before 1.0 Myr of cluster evolution. A break in the disc mass distributions is seen around 100 stars pc⁻², in particular for the R0.1 and R0.3 models, with masses dropping about one order of magnitude between 100 and 10⁴ stars pc⁻². Several surveys have shown that disc mass decreases in the vicinity of bright stars, and in regions with higher stellar number density (e.g. Fang et al. 2012; Mann & Williams 2009; Mann et al. 2014; Ansdell et al. 2017; van Terwisga et al. 2020). In Fig. 4, we show our simulation results and compare them with mean disc masses of various observed regions. It is important to note that, due to the nature of estimating disc masses from observations, the mean disc masses from our simulations and from observed regions differ systematically. In observations, mean disc masses are estimated by fitting a lognormal distribution to the measured dust masses (e.g. Williams et al. 2019). By calculating the same value for our simulations simply using the mean disc mass, the simulation curves are biased toward high-mass discs. Still, the behaviour observed in the simulated disc masses follows the trend of the observations: disc masses decrease as stellar density increases.

The trend in the disc mass distribution for local stellar densities 100–10⁴ stars pc⁻² suggests that, in our models, discs in regions in that density range are less likely to survive long enough or to have enough mass to form planetary systems. The planet formation process should already be underway before 1.0 Myr (Fig. 2) for discs to have the minimum reservoir of 10 M_⊕ in solids proposed by Ansdell et al. (2016) as necessary to form rocky planets or the cores of gas giants. From Fig. 4, it can be seen that, in our simulations, all discs in the R5.0, R2.5, and R1.0 models that are in areas of projected local density lower than 100 stars pc⁻² have masses in excess of 10 M_⊕. The mean disc dust mass in all other models is below this threshold by 2.0 Myr.

Most of the surviving discs in our simulations are around massive stars (0.5–1.9 M_⊙, see Fig. 2). A big factor in this is simply the construction of our models, where initial disc mass is proportional to stellar mass. Fig. 2, however, shows that drops of around 50 per cent in fractions of discs around high-mass stars are still present in high-density regions.

At the end of the simulations, the most massive discs are located in areas where the local stellar density is below 100 stars pc⁻³. This implies that large, massive discs observed today either formed in low-density regions or migrated to the outskirts of their birth locations fairly quickly. Discs born in the periphery of such regions have a much larger chance of surviving, and we could argue that the disc distributions seen in these low-density regions are similar to primordial disc distributions as they are pretty much unperturbed by external photoevaporation.

5 CONCLUSIONS

We perform simulations of star clusters with circumstellar discs. We implement the stellar dynamics, stellar evolution, viscous evolution of the discs, and external photoevaporation process to evolve simultaneously. We model our clusters as Plummer spheres with 10³ stars and initial virial radii of 0.1, 0.3, 0.5, 1.0, 2.5, and 5.0 pc to span a range of different number densities. Stars with masses M_* ≤ 1.9 M_⊙ are initially surrounded by a circumstellar disc, and stars with masses M_* > 1.9 M_⊙ do not have discs and are considered as only emitting UV radiation. Each cluster is evolved for 2.0 Myr. We can summarize our findings as follows:

1. External photoevaporation is efficient in destroying circumstellar discs quickly in all simulation models.

2. Mean disc lifetimes range from |${0.38 \pm 0.47}\, {\rm Myr}$| in the denser models (R_cluster = [0.1, 0.3, 0.5] pc), to |${0.65 \pm 0.55}\, {\rm Myr}$| for the sparser models (R_cluster = [1.0, 2.5, 5.0] pc).

3. Disc half-life, the time that it takes for half of the discs to be destroyed in a simulation run, ranges from |${0.20 \pm 0.01}\, {\rm Myr}$| in the denser models to |${1.42 \pm 0.33}\, {\rm Myr}$| in the sparser models.

4. Disc lifetimes, disc half-lives, disc fractions, and disc masses decrease as the stellar density of the models increase.

5. For the final disc masses in the denser regions (R_cluster = [0.1, 0.3, 0.5] pc), a projected local number density of 100 stars pc⁻² introduces a break in the distributions. There are only small variations in the masses of discs around stars in areas of lower densities. As the density increases beyond 100 stars pc⁻², the denser regions present a drop of almost an order of magnitude in disc masses.

6. The trends obtained in our simulations between disc mass and local stellar density are in agreement with dust mass measurements of discs in different observed regions: we compare our simulation results to masses of dusty young stellar objects in the Lupus clouds, the ONC, the Orion Molecular Cloud-2, Taurus, and NGC 2024.

ACKNOWLEDGEMENTS

We would like to thank the anonymous referee for their constructive comments, which greatly helped improve this paper. FC-R would like to thank the Star formation and protoplanetary disc group at Leiden Observatory for helpful discussions and insights. This work was carried out on the Dutch national e-infrastructure with the support of SURF Cooperative. This work was performed using resources provided by the Academic Leiden Interdisciplinary Cluster Environment (ALICE). This paper makes use of the packages numpy (Van Der Walt, Colbert & Varoquaux 2011), scipy (Virtanen et al. 2019), matplotlib (Hunter 2007), and makecite (Price-Whelan, Mechev & jumeroag 2018).

DATA AVAILABILITY

The data underlying this article are available at https://doi.org/10.5281/zenodo.3897171.

Footnotes

REFERENCES