External photoevaporation of circumstellar disks constrains the timescale for planet formation

Planet-forming circumstellar disks 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 disks are exposed to can vary in strength by at least six orders of magnitude. This radiation can limit the masses and sizes of the disks. Diversity in star forming environments can have long lasting effects in disk evolution and in the resulting planetary populations. We perform simulations to explore the evolution of circumstellar disks 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 disks: in regions of stellar density $\rho \sim 100 \mathrm{\ M}_\odot \mathrm{\ pc}^{-3}\mathrm{\ }$ around 80% of disks are destroyed before 2 Myr of cluster evolution. Our findings are in agreement with observed disk fractions in young star forming regions and support previous estimations that planet formation must start in timescales<0.1 - 1 Myr.


INTRODUCTION
Circumstellar disks 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 et al. 2010), during their first stages of evolution the disks remain embedded in an environment that is dense in gas and neighbouring stars. These conditions can be hostile for the disks in a myriad of ways: they can be subject to dynamical truncations (Vincke et al. 2015;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;, 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 disks through ram pressure stripping (Wijnen et al. 2016(Wijnen et al. , 2017. Since planet formation related processes seem to start very quickly in circumstellar disks (< 0.1 − 1 Myr, Najita & Kenyon (2014); Manara et al. (2018)), understanding the mechanisms that affect disk evolution is directly connected to understanding planetary E-mail: fconcha@strw.leidenuniv.nl 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 disks can help us comprehend the origins of the Solar System.
There are several observational indications that the environment of circumstellar disks shortly after their formation is unfavorable for their survival. Disks 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 et al. 2019). Moreover, observations indicate that disk fractions decline in regions close to an O-type star (e.g. Balog et al. 2007;Guarcello et al. 2007Guarcello et al. , 2009Guarcello et al. , 2010Fang et al. 2012;Mann et al. 2014;Guarcello et al. 2016). Fatuzzo & Adams (2008) estimate an FUV radiation field of up to G 0 ≈ 1000 in star clusters of N > 1000 stars 1 , while Facchini et al. (2016) show that disks of radius ∼ 150 au are subject to photoevaporation even in very low FUV fields (G 0 = 30). In regions of high stellar density, nearby stars can also affect disk size and morphology by dynamical interactions. Observational evidence of the imprints of dynamical truncations has been reported in several nearby stellar clusters (Olczak et al. 2 Concha-Ramírez et al. 2008;Reche et al. 2009;de Juan Ovelar et al. 2012). Tidal stripping that can be explained by disk-star interactions has been observed in the RW Aurigae system (Cabrit et al. 2006;Rodriguez et al. 2018;Dai et al. 2015) 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 disks 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 disks. As a consequence mass is accreted from the innermost regions of the disk onto 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 disks in the Lupus clouds versus disks 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 disks. External photoevaporation has been modeled with radiation hydrodynamics codes that solve flow equations through the disk boundaries, together with photodissociation region codes to obtain the temperature profiles of the disks (e.g. Haworth et al. 2016;Facchini et al. 2016). This method has been coupled with α-disk models to account for viscous evolution of the disk (e.g. Adams et al. 2004;Anderson et al. 2013;Gorti et al. 2015;Rosotti et al. 2017). Haworth et al. (2018a) introduce the concept of pre-computing photevaporation mass losses in terms of the surface density of the disks, an approach that we expand on in section 2.3. Winter et al. (2019a) model the environment of Cygnus OB2 and use the photoevaporation mass loss rate to constrain the timescale 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 disks (< 2 Myr) in moderate and low density regions ( 100 M pc −3 ).
Regarding dynamical effects, close encounters on a single Nbody disk of test particles have been investigated in several studies (e.g. Breslau et al. 2014;Jílková et al. 2016;Bhandare et al. 2016;Pfalzner et al. 2018). Winter et al. (2018a,b) 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 disk, and then use them to simulate the disk using a smoothed particles hydrodynamics (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 et al. 2006Olczak et al. , 2010Malmberg et al. 2011;Steinhausen & Pfalzner 2014;Vincke et al. 2015;Vincke & Pfalzner 2016. Directly adding SPH disks 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 M stars in a cluster with 100 stars. This allows them to study the effects of viscous spreading of the disks 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 et al. 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 disks inside stars clusters using a parametrized model for the disks. They show that viscous evolution and dynamical encounters are unable to explain the compact disks observed in star forming regions. They argue that other processes must affect the early evolution of the disks. Here we expand the model by improving the description of the viscous disks and by adding external photoevaporation due to the presence of bright nearby stars.
We model the circumstellar disks 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 disk. Dynamical truncations are parametrized, and the mass loss due to external photevaporation is calculated using the Far-ultraviolet Radiation Induced Evaporation of Discs (FRIED) grid ). This grid consists of pre-calculated mass loss rates for disks of different sizes and masses, immersed in several different external FUV fields. We use the Astrophysical Multipurpose Software Environment (AMUSE 2 , 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 3 .

Viscous growth of circumstellar disks
We implement circumstellar disks using the Viscous Accretion Disk Evolution Resource (VADER) by Krumholz & Forbes (2015).
VADER is an open source viscous evolution code which uses finitevolume discretization to solve the equations of mass transport, angular momentum, and internal energy in a thin, axisymmetric disk. An AMUSE interface for VADER was developed in the context of this work and is available online.
For the initial disk column density we use the standard disk profile introduced by Lynden- Bell & Pringle (1974): with where r c is the characteristic radius of the disk, r d and m d are the initial radius and mass of the disk, respectively, and Σ 0 is a normalization constant. Considering r c ≈ r d (Anderson et al. 2013), the density profile of the disk is: This expression allows the disk 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 disks we set the column density outside r d to a negligible value Σ edge = 10 −12 g cm −2 . The FRIED grid that we use to calculate the photoevaporation mass loss (see section 2.3.2) is a function of disk radius and outer surface density. There is a numerical challenge in determining what the disk 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 disk radius as the position of the first cell next to Σ edge , as shown in Figure 1.
The temperature profile of the disks is given by where T m is the midplane temperature at 1 au. Based on Anderson et al. (2013) we adopt T m = 300 K. Each disk is composed of a grid of 50 logarithmically spaced cells, in a range between 0.5 and 5000 au. In Appendix A we show that the resolution is enough for our calculations. The disks have a Keplerian rotation profile and turbulence parameter α = 5 × 10 −3 . The fact that the outer radius of the grid is much larger than the disk sizes (which were initially around 100 au, see section 2.4.1) allows the disks 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 disk radius. The mass flow through the inner boundary is considered as accreted mass and added to the mass of the host star.

Dynamical truncations
A close encounter between disks induces a discontinuity in their evolution. To modify the disks we calculate parametrized truncation radii. For two stars of the same mass, Rosotti et al. (2014) approximated the truncation radius to a third of the encounter distance. Together with the mass dependence from Bhandare et al. (2016), the dynamical truncation radius takes the form: where m 1 and m 2 are the masses of the encountering stars.
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 disk. To define r as the new disk radius we change the column density of all the disk cells outside r to the edge value Σ edge = 10 −12 g cm −2 .

External photoevaporation
A circumstellar disk can be evaporated by radiation coming from its host star or from a bright star in the vicinity. The heating of the disk surface can lead to the formation of gaps at different locations, which can cause the progression to a transition disk (e.g. Clarke et al. 2001;). The radiation can also truncate the disk 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 au of the disk (Font et al. 2004;Owen et al. 2010). Radiation from the host star can evaporate material from the outskirts of the disk, however Owen et al. (2010) show that more than 50% of the total mass loss occurs in the 5 − 20 au region. In Font et al. (2004) around 90% of the mass loss occurs within the inner 18 au of the disk. Given that the mass loss rate from the outer disk is typically comparable to or larger than that from the inner disk, we expect external photoevaporation to dominate in the outer regions. External photoevaporation has been shown to be more effective in evaporating the disks than radiation from the host star (e.g. Johnstone et al. 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 disks (Rosotti et al. 2017). In this work we ignore the effects of photoevaporation on internal disk structure, and deal exclusively with disk 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 extremeultraviolet (EUV), far-ultraviolet (FUV), and X-rays. In the case of external photoevaporation the dispersal of disk material is dominated by the FUV photons (Armitage 2000;Adams et al. 2004;). The main part of our work deals with photoevaporation due to FUV photons; in addition we also incorporate the effect of EUV photons (see Eq. 7).
The amount of mass lost from the disks 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 disks, is used to obtain the amount of radiation received by each disk. 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.

FUV luminosities
We follow Adams et al. (2004) in defining the FUV band ranging from 6 eV up 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 black body 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 precompute a relation between stellar mass and FUV luminosity. We do this by considering Figure 2. FUV luminosity vs. stellar mass fit calculated from the ZAMS spectra. M * = 1.9 M is the lower mass limit of the FRIED grid, and the lower mass limit for the stars to be considered emitting FUV radiation in our simulations (see Section 2.3.2 for details).
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 et al. 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 Figure 2. A similar fit was obtained by Parravano et al. (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 Figure 2 is 0.12 − 100 M , however, we are only interested in the range 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 disks 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 precomputed fit, after evolving the stellar evolution code. The low mass stars are not subject to stellar evolution.

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 disks. It spans a five dimension parameter space consisting of disk sizes (1 − 400 au), disk masses ( 10 −4 − 10 2 M Jup ), FUV fields (10 − 10 4 G 0 ), stellar masses (0.05 − 1.9 M ) and disk outer surface densities. The seemingly three dimensional grid subspace of disk mass, edge surface density, and disk radius is in fact two dimensional, as any combination of disk radius and disk 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, disk mass, and disk 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 disks 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 disk 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 disks do not always exactly match the ones in the grid. We perform interpolations over the grid to calculate the mass losses of the disks 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 disk in each dimension. The high resolution of the FRIED grid ensures that this interpolation is performed over an already smooth region.
When a massive star approaches a circumstellar disk, external photoevaporation is dominated by EUV radiation. Following Johnstone et al. (1998) we define a distance limit below which EUV photons dominate: where 2 fr Φ 49 1/2 ≈ 4, r d 14 = r d 10 14 cm with r d the disk radius, and 5 × 10 17 cm ∼ 3 × 10 4 au ∼ 0.16 pc. When a star with a disk is at a distance d < d min from a massive star we calculate the mass loss using equation (20) from Johnstone et al. (1998): with x ≈ 1.5 and ≈ 3. During most of their evolution, however, the circumstellar disks in the simulations experience photoevaporation only due to FUV photons. Since we do not consider interstellar gas and dust in the clusters, we do not account for extinction in the calculation of the radiation received by each small star.

Disk truncation due to photoevaporation
Once the mass loss due to photoevaporation is calculated for every disk, the disks 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 disk. We do this by moving outside-in from the disk edge and removing mass from each cell by turning its column density to the edge value Σ edge = 10 −12 defined in section 2.1. We stop at the point where the mass removed from the disk is equal to the calculated mass loss. We consider a disk to be completely evaporated when it has lost 99% of its initial mass (Anderson et al. 2013) or when its mean column density is lower than 1 g cm −2 (Pascucci et al. 2016). From this point forward the star continues its dynamical evolution normally, but is no longer affected by massive stars.

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 yr), internal evolutionary processes such as stellar evolution and gravitational dynamics proceed in their own internal time steps (500 yr and 1000 yr respectively). Throughout each macroscopic time step, we perform the following operations: 3) the code for the disk is stopped and removed and the star continues evolving only as part of the gravitational dynamics code. 6. Simulation runs until 5 Myr of evolution of until all the disks are dispersed, whichever happens first.
We present a scheme of this process and of the photoevaporation in Figures 3 and 4 respectively.

Disks
The initial radii of the circumstellar disks are given by: where R is a constant. The youngest circumstellar disks observed to date have diameters that range from ∼ 30 au (e.g. Lee et al. 2018) to ∼ 120 − 180 au (e.g. Murillo et al. 2013;van 't Hoff et al. 2018). Based on this we choose R = 100 au, which for our mass range 0.05 − 1.9 M for stars with disks results in initial disk radii between ∼ 22 au and ∼ 137 au. The initial masses of the disks are defined as:

Cluster
We perform simulations of young star clusters with stellar densities ρ ∼ 100 M pc −3 and ρ ∼ 50 M pc −3 using Plummer sphere spatial distributions (Plummer 1911). We will refer to these distributions as ρ 100 and ρ 50 , 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 M * ≈ 0.3 M . In Table 1 we show the details of the stellar masses in each simulation.
The simulations end at 5 Myr or when all the disks are dispersed, whichever happens first.

Disk 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 Figure  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 ρ 100 and ρ 50 regions, respectively. The shaded areas show the extent of the results in the different simulation runs. The mass lost from the circumstellar disks is dominated by external photoevaporation over the entire lifetime of the simulated clusters. Dynamical truncations only have a local effect on truncating disk radii and masses, whereas photoevaporation is a global effect influencing all disks in the cluster.
The amount of FUV radiation received by each disk 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 ρ 100 region the average FUV radiation over 5 Myr of evolution was ∼ 475 G 0 , with a minimum of 10 G 0 and a maximum of ∼ 10 4 G 0 . For the ρ 50 region the average over 5 Myr was ∼ 56 G 0 , minimum ∼ 2 G 0 and maximum 267 G 0 . These values are only shown as an indicative of the environment that the simulation disks 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 Figure 6 we show the cumulative distributions of disks 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 Figure 6 it can be seen that before 2 Myr of evolution the disks in ρ 100 lose mass more strongly than the ones in ρ 50 . However, starting around 2 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 ρ 100 the crossing time is t cross = 1.20 ± 0.04 Myr, and for ρ 50 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 Figure 7 we show the time step of maximum mass loss for each disk. 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 disk loses much mass. In Figure 8 we show how the cumulative distributions of disk masses at different moments in the simulation. The solid lines correspond to ρ 100 the dashed lines to ρ 50 . Each line corresponds to the total disks in all simulations. It can be seen that disk mass distributions decrease monotonically.
It is important to note that the FRIED grid has a lower limit of 10 G 0 for the FUV field, which is higher than the minimum experienced in the ρ 50 region. However, more than 90% 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 disk. We expand on the consequences of this for our results in section 4.3.

Disk lifetimes
The lifetime of circumstellar disks in young cluster regions is an important criterion to determine how photoevaporation affects planet formation. In Figure 9 we show disk fractions at different times of cluster evolution (black lines), together with observed disk 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.  The relaxation time is defined as:  (Spitzer 1987), where N is the number of stars, γ = 0.4, and the dynamical time is: where R is the radius of the cluster and M is its total mass. The relaxation time depends on the number of stars and the radius and mass of the stellar cluster. These are values that change through the dynamical evolution of a cluster, meaning that the relaxation time can grow and shrink at different time steps. These variations result in the jagged lines in Figure 9. For the simulations shown in Figure 9, t/t relax = 0.5 is reached at t = 2.01 ± 0.37 Myr of evolution for ρ 100 and at t = 2.05 ± 0.35 Myr for ρ 50 . Disk fractions in the ρ 100 simulations drop to around 20% before 2 Myr of cluster evolution. In the regions with ρ 50 the disks survive longer, but still most of the disks have disappeared by the end of the simulations.
Planet formation could still occur in disks that have been affected by photoevaporation as long as they are not completely dispersed. For gas giant cores and rocky planets to form, protoplanetary disks need to have a reservoir of dust mass M dust 10M ⊕ (Ansdell et al. 2018). In Figure 10 we show the fraction of disks with solid masses M disk > 10M ⊕ in time, for both simulated stellar densities. We use a 1:100 dust:gas mass ratio to turn the total disk mass into dust mass. For the ρ 100 regions the number of disks that fulfill this mass requirement drops to around 20% at 1 Myr, with less than 10% of disks of said mass still present at the end of the simulations. For the less dense regions, at the end of the simulations around 20% of disks with masses M disk > 10M ⊕ survive.
In order to make a parallel with the Solar System, in Figure 11 we show the number of disks in time with radii higher than 50 au, for both density regions. The drop in disk sizes is slower than the drop in disk masses as seen in Figure 10, and in the ρ 50 case the fraction of disks with radius > 50 au increases in the first time steps. This is related to the fact that, while some low mass disks get destroyed, others disks are still expanding due to viscous evolution. Some of these R disk > 50 au disks could still have masses or surface densities that are too low to form a planetary system.

Disk survival and consequences for planet formation
The results of the simulations carried out in this work characterize external photoevaporation as an important mechanism for disk dispersion. In comparison, the effect of dynamical truncations is negligible as a means for disk destruction.
The mean radiation received by the stars in our simulations fluctuates around ∼ 500 G 0 for the ρ 100 region and around ∼ 50 G 0 for the ρ 50 region. The FUV flux in the ONC is estimated to be ∼ 3 × 10 4 G 0 (O'dell & Wen 1994). Kim et al. (2016) estimate ∼ 3000 G 0 around a B star in NGC 1977, a region close to the ONC. According to this and to our results, most of the disks in such a dense region would be destroyed before reaching 1 Myr of age. Figure  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 0 and ∼ 4 G 0 , respectively. In particular, our results for the ρ 50 region show that even very low FUV fields can be effective in dispersing circumstellar disks over time. Winter et al. (2019b)  Protoplanetary disks 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. Figure 10 is in agreement with these conclusions. In our simulations, by 1 Myr around 20% of the disks have masses 10M ⊕ . This number drops to ∼ 10% by 3 Myr. According to our results rocky planets and gas giant cores must form very early on, otherwise the protoplanetary disks 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 disks are lower than the observed amount of heavy elements in extrasolar planetary systems around the same type of stars (see e.g. Manara et al. 2018;Najita & Kenyon 2014;Williams 2012;Greaves & Rice 2010, for discussions on this topic). Two scenarios have been proposed to explain this discrepancy in disk and exoplanet masses. The first one suggests that planet cores emerge within the first Myr of disk evolution, or even during the embedded phase while the disk is still being formed (e.g. Williams 2012; Greaves & Rice 2010). The second scenario proposes that disks can work as conveyor belts, transporting material from the surrounding interstellar medium towards the central star (e.g. Throop & Bally 2008;Kuffmeier et al. 2017).
Disk dispersal is not homogeneous across stellar types. There are observational indications that disk dispersion timescales depend on the mass of the host star, and that less massive (∼ 0.1 − 1.2 M ) stars keep their disks for longer than massive stars (Carpenter et al. 2006(Carpenter et al. , 2009Luhman & Mamajek 2012). We do not see this effect in our simulations, where the most massive stars keep their disks for longer simply because they initially have the most massive disks. The same effect is observed in Winter et al. (2019b), who used an analytic approach to estimate protoplanetary disk 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 disk dispersal. Radial drift of dust, fragmentation of large grains, and planetesimal formation are observed mechanisms that can affect both disk lifetimes and observed disk sizes. Viscous evolution alone is another internal process that can contribute to disk dispersal. A more complete model of disk evolution is needed to include the interplay between internal and external dispersion processes.
Initial disk masses are currently highly uncertain. Our chosen value of M d (t = 0) = 0.1M * is arbitrary, but disks of higher masses could still be stably supported (Haworth et al. 2018a;Nixon et al. 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 < 50% for planetary systems with eccentricities e ∼ 0.05 and semi-major axes < 20 au over 100 Myr of evolution. van Elteren et al. (2019) find that, in regions such as the Trapezium cluster, ∼ 30% of planetary systems are affected by the influence of other stars. Their fractal initial conditions provide local regions of higher densities, which are more favorable 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 disks, 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. Free-floating planets have been observed in the galactic center, and  efficiency analyses of these detections suggest that there are many more yet to be observed (e.g. Ryu et al. 2019;Navarro et al. 2019).

Influence of initial conditions
The effect of switching on photoevaporation when our simulations start have dramatic consequences for the initial circumstellar disks. Mass loss due to photoevaporation occurs very quickly once the stars are immersed in the FUV field. Around 60% and 20% of the disks are dispersed within the initial 50 000 yr in the ρ 100 and ρ 50 regions, respectively. The mean mass of the host stars whose disks dispersed in the initial 50 000 yr is 0.17 ± 0.03M for the ρ 100 region and 0.14 ± 0.04M for the ρ 50 region. In Figure 12 we show the disk fractions in time, separately for stars with masses M * < 0.5M and M * ≥ 0.5M . It can be seen that, for the stars of masses M * < 0.5M , the disk fractions drop much more dramatically during the first thousand years of cluster evolution.
In reality, if circumstellar disks 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 disks, increasing their lifetimes. They show that Cygnus OB2 probably underwent a primordial gas expulsion process that ended ∼ 0.5 Myr ago. This is based on the fact that 0.5 Myr of exposure to FUV fields reproduces the observed disk fractions in the region. Given that the estimated age of Cygnus OB2 is 3 − 5 Myr (Wright et al. 2010) the primordial gas in the region insulated the disks from external photoevaporation for several Myr. A similar point is made by Clarke (2007), who propose that the FUV field of the star θ 1 Orionis C must have been "switched on" no more than 1 − 2 Myr ago to explain the disk 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 disks and make the disk fraction drop more smoothly than what is shown in Figure 9.
The observations shown in Figure 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 disks. 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 disks 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 disk 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 disk 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 center of the cluster. This would increase the effect of photevaporation in the central regions of the cluster.

Model caveats
There are several physical processes not considered in this work which could affect the results presented here. One big caveat of our model is the lack of separation between dust and gas components in circumstellar disks. These separate disk components evolve differently and are affected in distinctive ways by outside mechanisms such as the ones implemented in this work. Gas disks has been observed to be larger than dust disks 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 disks 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 Haworth et al. 2018a) and consequences on planet formation. The conclusions regarding planet formation timescales derived in this work only consider the life expectancy of the disks, but considering different dust and gas disk components will likely affect these results.
While photoevaporation is considered to be primarily damaging for disks when coming from external sources, under certain regimes the photons coming from the host star can also contribute to disk dispersal.  and  show that FUV photons from the host star can drive photoevaporation mass loss at disk 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 ∼ 20 au region. Fatuzzo & Adams (2008) and Hollenbach et al. (2000) find that external photoevaporation dominates for disk regions > 10 au. The approach used in this work is valid for the disk 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 disk remained unperturbed. In reality, external FUV radiation can heat the whole surface of the disk, and mass loss can occur not just as a radial flow but also as a vertical flow from all over the disk (Adams et al. 2004). Given that the mass in the outer regions of a disk 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 1-dimensional disk model. New simulations by Haworth & Clarke (2019) show that, when considering a 2-dimensional disk model, mass loss rates can increase up to a factor 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 disks.
In the present 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 disks. Disks 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 disks are more compact and less bright than the ones around isolated stars. Disks around binary stars might also have shorter lifetimes, due to effects of the companion on the viscous timescale of the disk and also because of photoevaporation inside the system (Shadmehri et al. 2018;Rosotti & Clarke 2018).
Another process that can have important consequences in the evolution of circumstellar disks are supernovae explosions. Close & Pittard (2017) showed that nearby (0.3 pc) supernova explosions can cause mass loss rates of up to 1 × 10 −5 M yr −1 which can be sustained for about 200 yr. Only disks that are faced with the flow face-on manage to survive, but still lose 50% of their mass in the process. Portegies  show that a supernova explosion at a distance between 0.15 and 0.4 pc could create a misalignment of ∼ 5 • .6 between the star and its disk, which is consistent with the inclination of the plane of the Solar System. Such an event would also truncate a disk at around the edge of the Kuiper belt (42−55 au). Similar effects can be caused by other outcomes of stellar evolution, such as winds (Pelupessy & Portegies Zwart 2012).

CONCLUSIONS
We perform simulations of star clusters with stellar densities ρ ∼ 100 M pc −3 and ρ ∼ 50 M pc −3 . The stars with masses M * ≤ 1.9 M are surrounded by circumstellar disks. Stars with masses M * > 1.9 M are considered sufficiently massive stars to emit UV radiation, causing the disks around nearby stars to evaporate. The disks are subject to viscous growth, dynamical encounters, and external photoevaporation. The simulations span 5 Myr of cluster evolution. The main results of this work are: 1. In clusters with density ρ ∼ 100 M pc −3 around 80% of disks are destroyed by external photoevaporation within 2 Myr. The mean background FUV field is ∼ 500 G 0 . 2. In clusters with density ρ ∼ 50 M pc −3 around 50% of disks are destroyed by external photoevaporation within 2 Myr. The mean background FUV field is ∼ 50 G 0 . This shows that even very low FUV fields can be effective at destroying disks over long periods of time. 3. Mass loss caused by dynamical encounters is negligible compared to mass loss caused by external photoevaporation. Disk truncations that result from dynamical encounters are not an important process in setting observed disk size and mass distributions. 4. At 1 Myr, ∼ 20% of disks in the ρ ∼ 100 M pc −3 region and ∼ 50% of disks ρ ∼ 50 M pc −3 region have masses M disk ≥ 10 M ⊕ , the theoretical limit for gas giant core formation. 5. Our results support previous estimations that planet formation must start in timescales < 0.1 − 1 Myr (e.g. Najita & Kenyon 2014;Manara et al. 2018). 6. The obtained disk fractions in the different density regions, together with the quick dispersion of the disks in all the simulations, suggest that initial conditions are very important in the development of models of early protoplanetary disk evolution.  in this work were carried out in the Cartesius supercomputer, part of the Dutch national supercomputing facility SURFsara.

APPENDIX A: RESOLUTION OF THE DISKS
We use a resolution of 50 cells for the disks which gives us a good trade-off between computing time and acceptable results. In isolated disk evolution this causes an overestimation of disk radius by ∼ 10% on average over 1 Myr of disk evolution, compared with higher resolutions ( Figure A1). Since all the disks in the simulation are affected by photoevaporation from the start, no disks evolve as in the isolated case. Disk masses are overestimated by less than ∼ 5% compared to higher resolution runs ( Figure 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 disks. Given that we define a disk as dispersed when it has lost ∼ 90% of its initial mass, the slight mass overestimate obtained with the 50 cells resolution does not reflect in a quicker evaporation of the disks. This paper has been typeset from a T E X/L A T E X file prepared by the author.