External photoevaporation of circumstellar discs constrains the time-scale for planet formation

ABSTRACT

Planet-forming circumstellar discs are a fundamental part of the star formation process. Since stars form in a hierarchical fashion in groups of up to hundreds or thousands, the UV radiation environment that these discs are exposed to can vary in strength by at least six orders of magnitude. This radiation can limit the masses and sizes of the discs. Diversity in star forming environments can have long lasting effects in disc evolution and in the resulting planetary populations. We perform simulations to explore the evolution of circumstellar discs in young star clusters. We include viscous evolution, as well as the impact of dynamical encounters and external photoevaporation. We find that photoevaporation is an important process in destroying circumstellar discs: in regions of stellar density ρ ∼ 100 M_⊙ pc⁻³ around |$80{{\ \rm per\ cent}}$| of discs are destroyed before |$2\, \hbox{Myr}$| of cluster evolution. In regions of ρ ∼ 50 M_⊙ pc⁻³ around |$50{{\ \rm per\ cent}}$| of discs are destroyed in the same time-scale. Our findings are in agreement with observed disc fractions in young star-forming regions and support previous estimations that planet formation must start in time-scales <0.1–1 Myr.

1 INTRODUCTION

Circumstellar discs develop as a result of the star formation process (Williams & Cieza 2011). Since a non- negligible fraction of stars are not born in isolation (Bressert et al. 2010; King et al. 2012), and gas left over from the star formation process can linger for a few Myr (Portegies Zwart, McMillan & Gieles 2010), during their first stages of evolution the discs remain embedded in an environment that is dense in gas and neighbouring stars. These conditions can be hostile for the discs in a myriad of ways: they can be subject to dynamical truncations (Vincke, Breslau & Pfalzner 2015; Portegies Zwart 2016; Vincke & Pfalzner 2016) or be affected by processes related to stellar evolution, such as stellar winds (Pelupessy & Portegies Zwart 2012), supernovae explosions (Close & Pittard 2017; Portegies Zwart et al. 2018), and photoevaporation due to bright OB stars in the vicinity (e.g. Guarcello et al. 2016; Haworth et al. 2017). The surrounding gas can also shrink the discs through ram pressure stripping (Wijnen et al. 2016, 2017). Since planet formation related processes seem to start very quickly in circumstellar discs (⁠|$\lt 0.1\!-\!1\, \hbox{Myr}$|⁠; Najita & Kenyon 2014; Manara, Morbidelli & Guillot 2018), understanding the mechanisms that affect disc evolution is directly connected to understanding planetary system formation and survival. The Sun was probably born within a star cluster (Portegies Zwart et al. 2009), so discerning how the cluster environment affects the evolution of the discs can help us comprehend the origins of the Solar system.

There are several observational indications that the environment of circumstellar discs shortly after their formation is unfavourable for their survival. Discs have been observed to be evaporating in several young star-forming regions (e.g. Fang et al. 2012; de Juan Ovelar et al. 2012; Mann et al. 2014; Kim et al. 2016; van Terwisga, Hacar & van Dishoeck 2019). Moreover, observations indicate that disc fractions decline in regions close to an O-type star (e.g. Balog et al. 2007; Guarcello et al. 2007, 2009, 2010; Fang et al. 2012; Mann et al. 2014; Guarcello et al. 2016). Fatuzzo & Adams (2008) estimate an FUV radiation field of up to G₀ ≈ 1000 in star clusters of N > 1000 stars,¹ while Facchini, Clarke & Bisbas (2016) show that discs of radius |${\sim}150\, \hbox{au}$| are subject to photoevaporation even in very low FUV fields (G₀ = 30). In regions of high stellar density, nearby stars can also affect disc size and morphology by dynamical interactions. Observational evidence of the imprints of dynamical truncations has been reported in several nearby stellar clusters (Olczak, Pfalzner & Eckart 2008; Reche, Beust & Augereau 2009; de Juan Ovelar et al. 2012). Tidal stripping that can be explained by disc–star interactions has been observed in the RW Aurigae system (Cabrit et al. 2006; Dai et al. 2015; Rodriguez et al. 2018) and in the T Tauri binary system AS 205 (Salyk et al. 2014). There is also evidence that the Solar system might have been affected by a close encounter with another star during its early stages (Jílková et al. 2015; Pfalzner et al. 2018).

Circumstellar discs are not only affected by external processes but also by their internal viscous evolution. The combination of outwardly decreasing angular velocity together with outwardly increasing angular momentum causes shear stress forces inside the discs. As a consequence mass is accreted from the innermost regions of the disc on to its host star, whereas the outermost regions expand (Lynden-Bell & Pringle 1974). Tazzari et al. (2017) propose that the measured offsets in sizes and masses of discs in the Lupus clouds versus discs in the Taurus-Auriga and Ophiuchus regions can be explained as observational evidence of viscous spreading. However, Rosotti et al. (2019) argue that current surveys do not yet have the sufficient sensitivity to properly detect this phenomenon.

Different approaches have been implemented to study the effects of these processes on the lifetime of circumstellar discs. External photoevaporation has been modelled with radiation hydrodynamic codes that solve flow equations through the disc boundaries, together with photodissociation region codes to obtain the temperature profiles of the discs (e.g. Facchini et al. 2016; Haworth et al. 2016). This method has been coupled with α-disc models to account for viscous evolution of the disc (e.g. Adams et al. 2004; Anderson, Adams & Calvet 2013; Gorti, Hollenbach & Dullemond 2015; Rosotti et al. 2017). Haworth et al. (2018a) introduce the concept of pre-computing photoevaporation mass-losses in terms of the surface density of the discs, an approach that we expand on in Section 2.3. Winter, Clarke & Rosotti (2019a) model the environment of Cygnus OB2 and use the photoevaporation mass-loss rate to constrain the time-scale for gas expulsion in the young star-forming region. Nicholson et al. (2019) perform simulations of star-forming regions where FUV photoevaporation is implemented in post-processing, and find a very short lifetime for the discs (⁠|$\lt 2\, \hbox{Myr}$|⁠) in moderate and low-density regions (⁠|${\lesssim}{100}\, {\mathrm{M}_\odot \mathrm{\ pc}^{-3}}$|⁠).

Regarding dynamical effects, close encounters on a single N-body disc of test particles have been investigated in several studies (e.g. Breslau et al. 2014; Bhandare, Breslau & Pfalzner 2016; Jílková et al. 2016; Pfalzner et al. 2018). Winter et al. (2018a) and Winter et al. (2018b) use a ring of test particles around a star to obtain linearized expressions of the effect that a stellar encounter can have on the mass and morphology of the disc, and then use them to simulate the disc using a smoothed particle hydrodynamic (sph) code. A different approach for studying these effects is evolving the stellar dynamics of the cluster separately, and applying the effects of dynamical encounters afterwards (e.g. Olczak, Pfalzner & Spurzem 2006; Olczak, Pfalzner & Eckart 2010; Malmberg, Davies & Heggie 2011; Steinhausen & Pfalzner 2014; Vincke et al. 2015; Vincke & Pfalzner 2016, 2018). Directly adding sph discs to a simulation of a massive star cluster is computationally expensive, since a high resolution is needed over long time-scales. The closest effort corresponds to the work by Rosotti et al. (2014), in which individual sph codes are coupled to half of the |${1}\, {\mathrm{M}_\odot }$| stars in a cluster with 100 stars. This allows them to study the effects of viscous spreading of the discs and dynamical truncations in a self-consistent way, but they are limited by the computational resources needed for this problem. Parametrized approaches have also been developed, where the cluster dynamics and effects of truncations (Portegies Zwart 2016) and viscous spreading (Concha-Ramírez, Vaher & Portegies Zwart 2019) are considered simultaneously.

Concha-Ramírez et al. (2019) investigate the effect of viscous growth and dynamical truncations on the final sizes and masses of protoplanetary discs inside stars clusters using a parametrized model for the discs. They show that viscous evolution and dynamical encounters are unable to explain the compact discs observed in star-forming regions. They argue that other processes must affect the early evolution of the discs. Here, we expand the model by improving the description of the viscous discs and by adding external photoevaporation due to the presence of bright nearby stars.

We model the circumstellar discs using the Viscous Accretion Disk Evolution Resource (vader) (Krumholz & Forbes 2015). This code solves the equations of angular momentum and mass transport in a thin, axisymmetric disc. Dynamical truncations are parametrized, and the mass-loss due to external photoevaporation is calculated using the Far-ultraviolet Radiation Induced Evaporation of Discs (FRIED) grid (Haworth et al. 2018b). This grid consists of pre-calculated mass-loss rates for discs of different sizes and masses, immersed in several different external FUV fields. We use the Astrophysical Multipurpose Software Environment (amuse,² Portegies Zwart & McMillan 2018) framework to couple these codes along with cluster dynamics and stellar evolution. All the code developed for the simulations, data analyses, and figures of this paper is available in Github.³

2 MODEL

2.1 Viscous growth of circumstellar discs

We implement circumstellar discs using the Viscous Accretion Disk Evolution Resource (vader) by Krumholz & Forbes (2015). vader is an open source viscous evolution code that uses finite-volume discretization to solve the equations of mass transport, angular momentum, and internal energy in a thin, axisymmetric disc. An amuse interface for vader was developed in the context of this work and is available online.

This expression allows the disc to expand freely at the outer boundary while keeping the condition of zero torque at the inner boundary r_i.

To establish the radius of the discs we set the column density outside r_d to a negligible value |$\Sigma _{\mathrm{edge}} = 10^{-12}\, \mathrm{g\, cm}^{-2}$|⁠. The FRIED grid that we use to calculate the photoevaporation mass-loss (see Section 2.3.2) is a function of disc radius and outer surface density. There is a numerical challenge in determining what the disc outer surface density and radius actually are, since there is a large gradient in it down to the Σ_edge value. Computing a mass-loss rate for a very low outer surface density in this steep gradient would return an artificially low result. Considering this we define the disc radius as the position of the first cell next to Σ_edge, as shown in Fig. 1.

Each disc is composed of a grid of 50 logarithmically spaced cells, in a range between 0.5 and |$5000\, \hbox{au}$|⁠. In the Appendix, we show that the resolution is enough for our calculations. The discs have a Keplerian rotation profile and turbulence parameter α = 5 × 10⁻³. The fact that the outer radius of the grid is much larger than the disc sizes (which were initially around |$100\, \hbox{au}$|⁠, see Section 2.4.1) allows the discs to expand freely without reaching the boundaries of the grid. The mass flow through the outer boundary is set to zero in order to maintain the density Σ_edge needed to define the disc radius. The mass flow through the inner boundary is considered as accreted mass and added to the mass of the host star.

2.2 Dynamical truncations

To implement truncations, we first calculate the corresponding truncation radius caused by the encounter, according to equation (5). We consider all the mass outside this radius to be stripped from the disc. To define r′ as the new disc radius we change the column density of all the disc cells outside r′ to the edge value |$\Sigma _{\mathrm{edge}} = 10^{-12}\, \mathrm{g\, cm}^{-2}$|⁠.

2.3 External photoevaporation

A circumstellar disc can be evaporated by radiation coming from its host star or from a bright star in the vicinity. The heating of the disc surface can lead to the formation of gaps at different locations, which can cause the progression to a transition disc (e.g. Clarke, Gendrin & Sotomayor 2001; Gorti, Dullemond & Hollenbach 2009b). The radiation can also truncate the disc by removing material from the outer, loosely bound regions (e.g. Adams et al. 2004).

Models of internal photoevaporation have shown that most of the mass-loss in this case occurs in the inner |$20\, \hbox{au}$| of the disc (Font et al. 2004; Owen et al. 2010). Radiation from the host star can evaporate material from the outskirts of the disc; however, Owen et al. (2010) show that more than 50 per cent of the total mass-loss occurs in the |$5\!-\!20\, \hbox{au}$| region. In Font et al. (2004), around 90 per cent of the mass-loss occurs within the inner |$18\, \hbox{au}$| of the disc. Given that the mass-loss rate from the outer disc is typically comparable to or larger than that from the inner disc, we expect external photoevaporation to dominate in the outer regions. External photoevaporation has been shown to be more effective in evaporating the discs than radiation from the host star (e.g. Johnstone, Hollenbach & Bally 1998; Adams & Myers 2001). Truncation by external photoevaporation can also result in changes of the viscosity parameter α, which further affects the viscous evolution of the discs (Rosotti et al. 2017). In this work, we ignore the effects of photoevaporation on internal disc structure, and deal exclusively with disc survival rates. Because of this and the discussion above, we focus on external photoevaporation due to bright stars in the vicinity.

OB-type stars emit heating radiation in the form of extreme-ultraviolet (EUV), far-ultraviolet (FUV), and X-rays. In the case of external photoevaporation, the dispersal of disc material is dominated by the FUV photons (Armitage 2000; Adams et al. 2004; Gorti & Hollenbach 2009). The main part of our work deals with photoevaporation due to FUV photons; in addition, we also incorporate the effect of EUV photons (see equation 7).

The amount of mass lost from the discs as a result of external photoevaporation depends on the luminosity of the bright stars in the cluster. This luminosity, together with the distance from each of the massive stars to the discs, is used to obtain the amount of radiation received by each disc. We can then calculate the amount of mass lost. Below we explain what we consider to be massive stars and how we calculate the mass-loss rate.

2.3.1 FUV luminosities

We follow Adams et al. (2004) in defining the FUV band ranging from 6 to 13.6 eV, or approximately from 912 to 2070 Å. We calculate the FUV radiation of the stars in the simulations based on their spectral types. Given the presence of spectral lines in this band, we use synthetic stellar spectra rather than relying on blackbody approximations. The spectral library used is UVBLUE (Rodriguez-Merino et al. 2005), chosen for its high coverage of parameter space and high resolution, spanning the appropriate wavelength ranges. It spans a three-dimensional parameter space of stellar temperature, metallicity, and surface gravity.

We use the UVBLUE spectral library to pre-compute a relation between stellar mass and FUV luminosity. We do this by considering all the stars in the cluster to have solar metallicity (Z = 0.02). We then select the temperature and surface gravity spectra closer to the zero-age main-sequence value of each star, according to the parametrized stellar evolution code SeBa (Portegies Zwart & Verbunt 1996; Toonen, Nelemans & Zwart 2012). Given that the masses and radii of the stars are known, using the chosen spectra we build a relationship between stellar mass and FUV luminosity. This relation takes the shape of a segmented power law, as is shown in Fig. 2. A similar fit was obtained by Parravano, Hollenbach & McKee (2003). In runtime, the stars are subject to stellar evolution and the FUV luminosity for each star was calculated directly from this fit using the stellar mass.

The mass range of the fit in Fig. 2 is 0.12–100 M_⊙; however, we are only interested in the range of 1.9–50 M_⊙. As is further explained in Section 2.3.2, we only consider stars with masses higher that 1.9 M_⊙ to be emitting FUV radiation, and 50 M_⊙ is the theoretical upper limit for the stellar mass distribution. Stars with masses ≤1.9 M_⊙ are given discs and are affected by the massive stars. The massive stars are also subject to stellar evolution, which is implemented with SeBa through the amuse interface. The FUV luminosity of each massive star is calculated in every time-step, from the pre-computed fit, after evolving the stellar evolution code. The low-mass stars are not subject to stellar evolution.

2.3.2 Mass-loss rate due to external photoevaporation

To calculate the mass-loss due to FUV external photoevaporation we use the Far-ultraviolet Radiation Induced Evaporation of Discs (FRIED) grid developed by Haworth et al. (2018b). The FRIED grid is an open access grid of calculations of externally evaporating circumstellar discs. It spans a five-dimensional parameter space consisting of disc sizes (1–400 au), disc masses ( 10⁻⁴–10² M_Jup), FUV fields (10–10⁴ G₀), stellar masses (0.05–1.9 M_⊙), and disc outer surface densities. The seemingly three-dimensional grid subspace of disc mass, edge surface density, and disc radius is in fact two-dimensional, as any combination of disc radius and disc mass has only one edge density associated with it. Because of this, we only take into account a four-dimensional grid of radiation field strength, host star mass, disc mass, and disc radius.

Following the stellar mass ranges of the FRIED grid, we separate the stars in the simulations into two subgroups: massive stars and low-mass stars. Massive stars are all stars with initial masses higher than 1.9 M_⊙, while low-mass stars have masses ≤1.9 M_⊙. Only the low-mass stars have circumstellar discs in the simulations. The massive stars are considered as only generating FUV radiation and affecting the low-mass stars. In this way, we make sure that we stay within the stellar mass limits of the FRIED grid. Low-mass stars (≲1 M_⊙) have a negligible UV flux (Adams et al. 2006), so this approximation holds well for our purposes. Calculation of the FUV radiation emitted by the massive stars is further explained in Section 2.3.1. These star subgroups are considered only for the calculation of FUV radiation and photoevaporation mass-loss. In the gravity evolution code, the two subgroups are undistinguishable.

The FRIED grid allows to take a circumstellar disc with a specific mass, size, and density, around a star with a certain mass, and from the FUV radiation that it receives, obtain the photoevaporation mass-loss. However, the parameters of the simulated discs do not always exactly match the ones in the grid. We perform interpolations over the grid to calculate the mass-losses of the discs in the simulations. Because of computational constraints, we perform the interpolations on a subspace of the grid, such that it contains at least one data point above and below the phase-space point of the disc in each dimension. The high resolution of the FRIED grid ensures that this interpolation is performed over an already smooth region.

2.3.3 Disc truncation due to photoevaporation

Once the mass-loss due to photoevaporation is calculated for every disc, the discs are truncated at a point that coincides with the amount of mass lost in the process. We take the approach of Clarke (2007) and remove mass from the outer edge of the disc. We do this by moving outside-in from the disc edge and removing mass from each cell by turning its column density to the edge value Σ_edge = 10⁻¹² defined in Section 2.1. We stop at the point where the mass removed from the disc is equal to the calculated mass-loss.

We consider a disc to be completely evaporated when it has lost 99 per cent of its initial mass (Anderson et al. 2013) or when its mean column density is lower than 1 g cm⁻² (Pascucci et al. 2016). From this point forward the star continues its dynamical evolution normally, but is no longer affected by massive stars.

2.3.4 Summary of cluster evolution

The code runs in major time-steps, which represent the time-scale on which the various processes are coupled. Within each of these macroscopic time-steps (⁠|$1000\, \hbox{yr}$|⁠), internal evolutionary processes such as stellar evolution and gravitational dynamics proceed in their own internal time-steps (⁠|$500\, \hbox{yr}$| and |$1000\, \hbox{yr}$|⁠, respectively). Throughout each macroscopic time-step, we perform the following operations:

1. Gravitational dynamics code is evolved.

2. We check the stars for dynamical encounters:

   • If a dynamical encounter occurs, we determine the truncation radius for each disc.

   • We update the radius and mass of the discs.

3. Stellar evolution code is evolved.

4. Photoevaporation process is carried out as follows. For each massive star,

   • We calculate the distance d from the massive star to each low-mass star.

   • If d < d_min (see equation 7), we calculate the mass-loss |$\dot{M}_{\rm EUV}$|⁠.

   • If d ≥ d_min the massive star’s FUV luminosity L_FUV is calculated (see Section 2.3.1).

   • Using d and L_FUV, we calculate l_FUV, the amount of FUV radiation received by the low-mass star.

   • Using l_FUV together with the low-mass star’s mass, disc mass, and disc radius, we build a subgrid of the FRIED grid and interpolate over it to calculate the mass-loss |$\dot{M}_{\rm FUV}$|⁠.

   • The total mass-loss in the time-step is calculated using |$\dot{M}_{\rm EUV}$| and |$\dot{M}_{\rm FUV}$|⁠.

   • The mass is removed from the disc by moving outside-in and removing mass from the cells.

   • The disc mass and radius are updated.

5. Discs are checked for dispersal. If a disc has been dispersed (see Section 2.3.3) the code for the disc is stopped and removed and the star continues evolving only as part of the gravitational dynamics code.

6. Simulation runs until |$5\, \hbox{Myr}$| of evolution of until all the discs are dispersed, whichever happens first.

We present a scheme of this process and of the photoevaporation in Figs 3 and 4, respectively.

2.4 Initial conditions

2.4.1 Discs

2.4.2 Cluster

We perform simulations of young star clusters with stellar densities ρ ∼ 100 M_⊙ pc⁻³ and ρ ∼ 50 M_⊙ pc⁻³ using Plummer sphere spatial distributions (Plummer 1911). We will refer to these distributions as ρ₁₀₀ and ρ₅₀ , respectively. All the regions are in virial equilibrium (viral ratio Q = 0.5).

Stellar masses are randomly drawn from a Kroupa mass distribution (Kroupa 2001) with maximum mass 50 M_⊙. The mean mass of the distribution is |$\overline{{M}}_* \approx 0.3 \mathrm{\ M}_\odot$|⁠. In Table 1, we show the details of the stellar masses in each simulation.

The simulations end at |$5\, \hbox{Myr}$| or when all the discs are dispersed, whichever happens first.

3 RESULTS

3.1 Disc mass-loss in time

As a way to quantify the mass-loss effect of each of the processes included in the simulations, we measure the mass-loss due to external photoevaporation and dynamical truncations separately. In Fig. 5, we show the mass-loss over time for external photoevaporation and dynamical truncations. The solid and dashed lines correspond to the mean mass-loss among all stars in the ρ₁₀₀ and ρ₅₀ regions, respectively. The shaded areas show the extent of the results in the different simulation runs.

The mass lost from the circumstellar discs is dominated by external photoevaporation over the entire lifetime of the simulated clusters. Dynamical truncations only have a local effect on truncating disc radii and masses, whereas photoevaporation is a global effect influencing all discs in the cluster.

The amount of FUV radiation received by each disc and the ensuing mass-loss are variable. The effect depends on the proximity to massive stars, which changes with time as the stars orbit in the cluster potential. For the ρ₁₀₀ region, the average FUV radiation over |$5\, \hbox{Myr}$| of evolution was ∼475 G₀ , with a minimum of 10 G₀ and a maximum of ∼10⁴ G₀ . For the ρ₅₀ region, the average over |$5\, \hbox{Myr}$| was ∼56 G₀ , minimum ∼2 G₀, and maximum 267 G₀ . These values are only shown as an indicative of the environment that the simulation discs were dispersed in; however, a short exposure to a strong FUV field can be instantly more destructive than a sustained low FUV field. The FUV field can also vary in time due to processes intrinsic to star formation, such as a massive star being embedded during its early stages (e.g. Ali & Harries 2019). In our simulations, however, photoevaporation is driven by the background FUV field. In Fig. 6 we show the cumulative distributions of discs that lose more than 5 M_Jup in time. It can be seen that large mass-losses are not driven by close encounters with bright stars but by the environmental FUV radiation. From Fig. 6 it can be seen that, before |$2\, \hbox{Myr}$| of evolution, the discs in ρ₁₀₀ lose mass more strongly than the ones in ρ₅₀ . However, starting around |$2\, \hbox{Myr}$| and until the end of the simulations there is large scatter in the mass-loss behaviour for each region. This is related to the dynamics of each cluster. For ρ₁₀₀ the crossing time is t_cross = 1.20 ± 0.04 Myr, and for ρ₅₀, the crossing time is t_cross = 0.98 ± 0.09 Myr. This results in the fact that, after one crossing time, all the stars in the clusters have been affected by the radiating stars similarly, causing the scatter in the mass-loss effects. While the density of each cluster defines the effects of photoevaporation in the early evolutionary stages, after one crossing time the initial density of the region is not as important and photoevaporation works uniformly.

In Fig. 7, we show the time-step of maximum mass-loss for each disc. It can be seen that, other than the effect of switching on the simulation at the beginning (see Section 4.2), there is not a particular time at which a disc loses much mass. In Fig. 8, we show how the cumulative distributions of disc masses at different moments in the simulation. The solid lines correspond to ρ₁₀₀ the dashed lines to ρ₅₀ . Each line corresponds to the total discs in all simulations. It can be seen that disc mass distributions decrease monotonically.

It is important to note that the FRIED grid has a lower limit of 10 G₀ for the FUV field, which is higher than the minimum experienced in the ρ₅₀ region. However, more than 90 per cent of the stars in these simulations are within the limits of the grid at all times. In the few cases where stars were outside the limits of the grid, the mass-loss obtained reflects a lower bound defined by the grid, but this does not affect our results.

Photoevaporation mass-loss can have different effects over the gas and dust components of a circumstellar disc. We expand on the consequences of this for our results in Section 4.3.

3.2 Disc lifetimes

The lifetime of circumstellar discs in young cluster regions is an important criterion to determine how photoevaporation affects planet formation. In Fig. 9, we show disc fractions at different times of cluster evolution (black lines), together with observed disc fractions from Ribas et al. (2014) and Richert et al. (2018). The orange line shows the mean of the observations, calculated using a moving bin spanning 10 observation points. The calculation of the mean starts by binning the first 10 points, and then sliding horizontally through the observations one point at a time such that 10 points are always considered.

For the simulations shown in Fig. 9, t/t_relax = 0.5 is reached at |${t} = 2.01 \pm 0.37\, \hbox{Myr}$| of evolution for ρ₁₀₀ and at |${t} = 2.05 \pm 0.35\, \hbox{Myr}$| for ρ₅₀ . Disc fractions in the ρ₁₀₀ simulations drop to around |$20{{\ \rm per\ cent}}$| before |$2\, \hbox{Myr}$| of cluster evolution. In the regions with ρ₅₀ the discs survive longer, but still most of the discs have disappeared by the end of the simulations.

Planet formation could still occur in discs that have been affected by photoevaporation as long as they are not completely dispersed. For gas giant cores and rocky planets to form, protoplanetary discs need to have a reservoir of dust mass M_dust ≳ 10M_⊕ (Ansdell et al. 2018). In Fig. 10 we show the fraction of discs with solid masses M_disc > 10M_⊕ in time, for both simulated stellar densities. We use a 1:100 dust:gas mass ratio to turn the total disc mass into dust mass. For the ρ₁₀₀ regions, the number of discs that fulfil this mass requirement drops to around |$20{{\ \rm per\ cent}}$| at |$1\, \hbox{Myr}$|⁠, with less than |$10{{\ \rm per\ cent}}$| of discs of said mass still present at the end of the simulations. For the less dense regions, at the end of the simulations around |$20{{\ \rm per\ cent}}$| of discs with masses M_disc > 10M_⊕ survive.

In order to make a parallel with the Solar system (Fig. 11), we show the number of discs in time with radii higher than |$50\, \hbox{au}$| for both density regions. The drop in disc sizes is slower than the drop in disc masses as seen in Fig. 10, and in the ρ₅₀ case, the fraction of discs with radius |$\gt 50\, \hbox{au}$| increases in the first time-steps. This is related to the fact that, while some low-mass discs get destroyed, others discs are still expanding due to viscous evolution. Some of these |${R}_\mathrm{disc} \gt 50\, \hbox{au}$| discs could still have masses or surface densities that are too low to form a planetary system.

4 DISCUSSION

4.1 Disc survival and consequences for planet formation

The results of the simulations carried out in this work characterize external photoevaporation as an important mechanism for disc dispersion. In comparison, the effect of dynamical truncations is negligible as a means for disc destruction.

The mean radiation received by the stars in our simulations fluctuates around ∼500 G₀ for the ρ₁₀₀ region and around ∼50 G₀ for the ρ₅₀ region. The FUV flux in the ONC is estimated to be ∼3 × 10⁴ G₀ (O’dell & Wen 1994). Kim et al. (2016) estimate ∼3000 G₀ around a B star in NGC 1977, a region close to the ONC. According to this and to our results, most of the discs in such a dense region would be destroyed before reaching |$1\, \hbox{Myr}$| of age. Fig. 9 also agrees with results by Facchini et al. (2016) and Haworth et al. (2017), which show that photoevaporation mass-loss can be important even for regions with ∼30 G₀ and ∼4 G₀, respectively. In particular, our results for the ρ₅₀ region show that even very low FUV fields can be effective in dispersing circumstellar discs over time. Winter et al. (2019b) find similar dispersion time-scales, with a median of |$2.3\, \hbox{Myr}$| in the solar neighbourhood and |$0.5\, \hbox{Myr}$| in the central regions of the Milky Way, for stars down to |$0.8\, \mathrm{M}_\odot$|⁠. Comparable results are obtained by Nicholson et al. (2019), who find the half-life of protoplanetary discs to be around |$2\!-\!3\, \hbox{Myr}$| in clusters of various initial conditions.

Protoplanetary discs need to have a reservoir of dust mass M_dust ≳ 10M_⊕ to be able to form the rocky cores of giant gas planets (Ansdell et al. 2018). Manara et al. (2018) show that such cores need to be already in place at ages ∼1–3 Myr for this type of planets to form. Fig. 10 is in agreement with these conclusions. In our simulations, by |$1\, \hbox{Myr}$| around |$20{{\ \rm per\ cent}}$| of the discs have masses ≳10M_⊕. This number drops to |${\sim}10{{\ \rm per\ cent}}$| by |$3\, \hbox{Myr}$|⁠. According to our results rocky planets and gas giant cores must form very early on, otherwise the protoplanetary discs are not massive enough to provide the necessary amount of solids. This is in agreement with observational time constrains for planet formation and with the so-called missing-mass problem: solids mass measurements in protoplanetary discs are lower than the observed amount of heavy elements in extrasolar planetary systems around the same type of stars (see e.g. Greaves & Rice 2010; Williams 2012; Najita & Kenyon 2014; Manara et al. 2018; for discussions on this topic). Two scenarios have been proposed to explain this discrepancy in disc and exoplanet masses. The first one suggests that planet cores emerge within the first Myr of disc evolution, or even during the embedded phase while the disc is still being formed (e.g. Greaves & Rice 2010; Williams 2012). The second scenario proposes that discs can work as conveyor belts, transporting material from the surrounding interstellar medium towards the central star (e.g. Throop & Bally 2008; Kuffmeier, Haugbøle & Nordlund 2017).

Disc dispersal is not homogeneous across stellar types. There are observational indications that disc dispersion time-scales depend on the mass of the host star, and that less massive (∼0.1–1.2 M_⊙) stars keep their discs for longer than massive stars (Carpenter et al. 2006, 2009; Luhman & Mamajek 2012). We do not see this effect in our simulations, where the most massive stars keep their discs for longer simply because they initially have the most massive discs. The same effect is observed in Winter et al. (2019b), who used an analytic approach to estimate protoplanetary disc dispersal time-scales by external photoevaporation. The discrepancy between observations and theoretical results suggests that internal processes not considered in this work can also play an important role in disc dispersal. Radial drift of dust, fragmentation of large grains, and planetesimal formation are observed mechanisms that can affect both disc lifetimes and observed disc sizes. Viscous evolution alone is another internal process that can contribute to disc dispersal. A more complete model of disc evolution is needed to include the interplay between internal and external dispersion processes.

Initial disc masses are currently highly uncertain. Our chosen value of M_d(t = 0) = 0.1M_* is arbitrary, but discs of higher masses could still be stably supported (Haworth et al. 2018a; Nixon, King & Pringle 2018).

Once a planetary system has formed, its survival inside a star cluster is not guaranteed. Of the 4071 exoplanets confirmed to date, only 30 have been found inside star clusters. Cai et al. (2019) performed simulations of planetary systems in dense, young massive star clusters. They found that the survival rate is |$\lt 50{{\ \rm per\ cent}}$| for planetary systems with eccentricities e ∼ 0.05 and semimajor axes |$\lt 20\, \hbox{au}$| over |$100\, \hbox{Myr}$| of evolution. van Elteren et al. (2019) find that, in regions such as the Trapezium cluster, |${\sim}30{{\ \rm per\ cent}}$| of planetary systems are affected by the influence of other stars. Their fractal initial conditions provide local regions of higher densities, which are more favourable for dynamical encounters than our initial conditions. When making parallels with currently observed exoplanet systems, it is important not only to consider the environment effects on the early protoplanetary discs but also on the planets themselves once they are already formed.

Observations suggest that planets are able to circumvent all of these adversary processes and still form in highly unlikely regions. Evidence of star formation and even proplyd-like objects have been observed around Sgr A* (Yusef-Zadeh et al. 2015, 2017). Free-floating planets have been observed in the galactic centre, and efficiency analyses of these detections suggest that there are many more yet to be observed (e.g. Navarro et al. 2019; Ryu et al. 2019).

4.2 Influence of initial conditions

The effect of switching on photoevaporation when our simulations start have dramatic consequences for the initial circumstellar discs. Mass-loss due to photoevaporation occurs very quickly once the stars are immersed in the FUV field. Around |$60{{\ \rm per\ cent}}$| and |$20{{\ \rm per\ cent}}$| of the discs are dispersed within the initial |$50\, 000\, \hbox{yr}$| in the ρ₁₀₀ and ρ₅₀ regions, respectively. The mean mass of the host stars whose discs dispersed in the initial |$50\, 000\, \hbox{yr}$| is |$0.17 \pm 0.03\, \mathrm{M}_\odot$| for the ρ₁₀₀ region and |$0.14 \pm 0.04\, \mathrm{M}_\odot$| for the ρ₅₀ region. In Fig. 12, we show the disc fractions in time, separately for stars with masses |${\it M}_* \lt 0.5\, \mathrm{M}_\odot$| and |${\it M}_* \ge 0.5\, {\rm M}_\odot$|⁠. It can be seen that, for the stars of masses |${\it M}_* \,\lt 0.5\, \mathrm{M}_\odot$|⁠, the disc fractions drop much more dramatically during the first thousand years of cluster evolution.

In reality, if circumstellar discs do form around such low-mass stars, they could be sheltered from photoevaporation by interstellar gas and dust, which can linger for several million years after star formation (Portegies Zwart et al. 2010). Models of the Cygnus OB2 region by Winter et al. (2019a) demonstrate that the extinction of FUV photons through the gas dampens the mass-loss of the discs, increasing their lifetimes. They show that Cygnus OB2 probably underwent a primordial gas expulsion process that ended |${\sim}0.5\, \hbox{Myr}$| ago. This is based on the fact that |$0.5\, \hbox{Myr}$| of exposure to FUV fields reproduces the observed disc fractions in the region. Given that the estimated age of Cygnus OB2 is |$3\!-\!5\, \hbox{Myr}$| (Wright et al. 2010) the primordial gas in the region insulated the discs from external photoevaporation for several Myr. A similar point is made by Clarke (2007), who propose that the FUV field of the star θ₁ Orionis C must have been ‘switched on’ no more than |$1\!-\!2\, \hbox{Myr}$| ago to explain the disc-size distribution observed around it. This switching on could have been caused by the star clearing out the primordial gas it was embedded in, thus reducing extinction around it and making its effective radiation field stronger (Ali & Harries 2019). The presence of gas in young star clusters could then protect the protoplanetary discs and make the disc fraction drop more smoothly than what is shown in Fig. 9.

The observations shown in Fig. 9 span clusters of many different ages and densities. Our simulation results show that one order of magnitude difference in initial cluster density can yield an important difference in the number of surviving discs. A one order of magnitude spread in cluster density translates to a three order of magnitude difference in cluster radius. The extent of stellar densities in regions where circumstellar discs have been detected does not only reflect the environment of these regions, but also the variety in initial cluster densities.

The initial spatial distribution of the stars in the simulations also plays an important role during the early stages of disc evolution. The stars in our simulations were initially distributed in a Plummer sphere with a specified radius and in virial equilibrium. An approach with fractal or filamentary (e.g. Winter et al. 2019a) initial conditions could change the overall disc survival rates. If a massive star is born in a clump of a fractal distribution, for example, stars in other clumps without massive stars could be minorly affected by radiation and have higher chances of surviving and, eventually, making planets. Higher density regions also increase the relevance of dynamical truncations. This effect of initial conditions could also be counteracted by dynamical mass segregation, in which massive stars move towards the centre of the cluster. This would increase the effect of photoevaporation in the central regions of the cluster.

4.3 Model caveats

There are several physical processes not considered in this work that could affect the results presented here. One big caveat of our model is the lack of separation between dust and gas components in circumstellar discs. These separate disc components evolve differently and are affected in distinctive ways by outside mechanisms such as the ones implemented in this work. Gas discs has been observed to be larger than dust discs by a factor of ∼2 (Ansdell et al. 2018). Whether this is caused by different evolution for gas and dust or observational optical depth effects is still up for debate (see e.g. Birnstiel & Andrews 2013; Facchini et al. 2017; Trapman et al. 2019, for discussions on the topic). The dust in protoplanetary discs is subject to radial drifting and radially dependent grain growth, which can make it resilient to photoevaporation. This can have direct implications on the photoevaporation mass-loss rates (Facchini et al. 2016; Haworth et al. 2018a) and consequences on planet formation. The conclusions regarding planet formation time-scales derived in this work only consider the life expectancy of the discs, but considering different dust and gas disc components will likely affect these results.

While photoevaporation is considered to be primarily damaging for discs when coming from external sources, under certain regimes the photons coming from the host star can also contribute to disc dispersal. Gorti & Hollenbach (2009) and Gorti, Dullemond & Hollenbach (2009a) show that FUV photons from the host star can drive photoevaporation mass-loss at disc radius ∼1–10 au and ≳ 30 au. Owen et al. (2010) and Font et al. (2004) show that internal photoevaporation can also remove loosely bound material from the outer regions; however, the largest mass-loss was from the inner |${\sim}20\, \hbox{au}$| region. Fatuzzo & Adams (2008) and Hollenbach, Yorke & Johnstone (2000) find that external photoevaporation dominates for disc regions |$\gt 10\, \hbox{au}$|⁠. The approach used in this work is valid for the disc truncation approximation; however, a more complete analysis would have to consider the combined action and interplay of external and internal photoevaporation.

Mass-loss due to photoevaporation was modelled by calculating a truncation radius and removing all the mass outside it, while the inner region of the disc remained unperturbed. In reality, external FUV radiation can heat the whole surface of the disc, and mass-loss can occur not just as a radial flow but also as a vertical flow from all over the disc (Adams et al. 2004). Given that the mass in the outer regions of a disc is more loosely bound to its host star, the truncation approach is a good first-order approximation for mass-loss. Furthermore, the FRIED grid used to estimate the photoevaporation mass-loss was built using a one-dimensional disc model. New simulations by Haworth & Clarke (2019) show that, when considering a two-dimensional disc model, mass-loss rates can increase up to a factor of 3.7 for a solar mass star. The photoevaporation mass-losses obtained in this work should then be considered as lower limits, but are still a good estimate of the effects of bright stars in the vicinity of circumstellar discs.

In this work, we did not include binary stars or any multiples. The presence of multiple stellar systems can have direct consequences on the dynamical evolution of the cluster and on the effects of photoevaporation over the discs. Discs around binary stars have been observed in the star-forming regions ρ Ophiuchus (Cox et al. 2017) and Taurus-Auriga (Harris et al. 2012; Akeson & Jensen 2014; Akeson et al. 2019). Observations suggest that these discs are more compact and less bright than the ones around isolated stars. Discs around binary stars might also have shorter lifetimes, due to effects of the companion on the viscous time-scale of the disc and also because of photoevaporation inside the system (Rosotti & Clarke 2018; Shadmehri, Ghoreyshi & Alipour 2018).

Another process that can have important consequences in the evolution of circumstellar discs are supernovae explosions. Close & Pittard (2017) showed that nearby (⁠|$0.3\, \hbox{pc}$|⁠) supernova explosions can cause mass-loss rates of up to |$1 \times 10^{-5}\, \mathrm{M_\odot \, yr^{-1}}$| that can be sustained for about |$200\, \hbox{yr}$|⁠. Only discs that are faced with the flow face-on manage to survive, but still lose |$50{{\ \rm per\ cent}}$| of their mass in the process. Portegies Zwart et al. (2018) show that a supernova explosion at a distance between 0.15 and |$0.4\, \hbox{pc}$| could create a misalignment of ∼5.6° between the star and its disc, which is consistent with the inclination of the plane of the Solar system. Such an event would also truncate a disc at around the edge of the Kuiper belt (⁠|$42\!-\!55\, \hbox{au}$|⁠). Similar effects can be caused by other outcomes of stellar evolution, such as winds (Pelupessy & Portegies Zwart 2012).

5 CONCLUSIONS

We perform simulations of star clusters with stellar densities ρ ∼ 100 M_⊙ pc⁻³ and ρ ∼ 50 M_⊙ pc⁻³ . The stars with masses M_* ≤ 1.9 M_⊙ are surrounded by circumstellar discs. Stars with masses M_* > 1.9 M_⊙ are considered sufficiently massive stars to emit UV radiation, causing the discs around nearby stars to evaporate. The discs are subject to viscous growth, dynamical encounters, and external photoevaporation. The simulations span |$5\, \hbox{Myr}$| of cluster evolution. The main results of this work are as follows:

1. In clusters with density ρ ∼ 100 M_⊙ pc⁻³ around |$80{{\ \rm per\ cent}}$| of discs are destroyed by external photoevaporation within |$2\, \hbox{Myr}$|⁠. The mean background FUV field is ∼500 G₀.

2. In clusters with density ρ ∼ 50 M_⊙ pc⁻³ around |$50{{\ \rm per\ cent}}$| of discs are destroyed by external photoevaporation within |$2\, \hbox{Myr}$|⁠. The mean background FUV field is ∼50 G₀. This shows that even very low FUV fields can be effective at destroying discs over long periods of time.

3. Mass-loss caused by dynamical encounters is negligible compared to mass-loss caused by external photoevaporation. Disc truncations that result from dynamical encounters are not an important process in setting observed disc size and mass distributions.

4. At |$1\, \hbox{Myr}$|⁠, |${\sim}20{{\ \rm per\ cent}}$| of discs in the ρ ∼ 100 M_⊙ pc⁻³ region and |${\sim}50{{\ \rm per\ cent}}$| of discs ρ ∼ 50 M_⊙ pc⁻³ region have masses M_disc ≥ 10 M_⊕, the theoretical limit for gas giant core formation.

5. Our results support previous estimations that planet formation must start in time-scales <0.1–1 Myr (e.g. Najita & Kenyon 2014; Manara et al. 2018).

6. The obtained disc fractions in the different density regions, together with the quick dispersion of the discs in all the simulations, suggest that initial conditions are very important in the development of models of early protoplanetary disc evolution.

ACKNOWLEDGEMENTS

We would like to thank the anonymous referee for their thoughtful comments that helped improve this paper. We would also like to thank Andrew Winter, Sierk van Terwisga, and the protoplanetary disc group at Leiden Observatory for helpful discussions and comments. FC-R would like to thank Valeriu Codreanu from SURFsara for his invaluable technical assistance. The simulations performed in this work were carried out in the Cartesius supercomputer, part of the Dutch national supercomputing facility SURFsara.

Footnotes

REFERENCES

APPENDIX: RESOLUTION OF THE DISCS

We use a resolution of 50 cells for the discs that gives us a good trade-off between computing time and acceptable results. In isolated disc evolution, this causes an overestimation of disc radius by |${\sim}10{{\ \rm per\ cent}}$| on average over 1 Myr of disc evolution, compared with higher resolutions (Fig. A1). Since all the discs in the simulation are affected by photoevaporation from the start, no discs evolve as in the isolated case. Disc masses are overestimated by less than |${\sim}5{{\ \rm per\ cent}}$| compared to higher resolution runs (Fig. A2). This results, in turn, in a slight underestimation of the effects that mass removal, whether through photoevaporation or dynamical encounters, has on the survival times of the discs. Given that we define a disc as dispersed when it has lost |${\sim}90{{\ \rm per\ cent}}$| of its initial mass, the slight mass overestimate obtained with the 50 cells resolution does not reflect in a quicker evaporation of the discs.