The external photoevaporation of structured protoplanetary disks

Context. The dust in planet-forming disks is known to evolve rapidly through growth and radial drift. In the high irradiation environments of massive star-forming regions where most stars form, external photoevaporation also contributes to rapid dispersal of disks. This raises the question of why we still observe quite high disk dust masses in massive star-forming regions. Aims. We test whether the presence of substructures is enough to explain the survival of the dust component and observed millimeter continuum emission in protoplanetary disks located within massive star-forming regions. We also characterize the dust content removed by the photoevaporative winds. Methods. We performed hydrodynamical simulations (including gas and dust evolution) of protoplanetary disks subject to irradiation fields of F UV = 10 2 , 10 3 , and 10 4 G 0 , with different dust trap configurations. We used the FRIED grid to derive the mass loss rate for each irradiation field and disk properties, and then proceed to measure the evolution of the dust mass over time. For each simulation we estimated the continuum emission at λ = 1 . 3mm along with the radii encompassing 90% of the continuum flux, and characterized the dust size distribution entrained in the photoevaporative winds, in addition to the resulting far-ultraviolet (FUV) cross section. Results. Our simulations show that the presence of dust traps can extend the lifetime of the dust component of the disk to a few millionyears if the FUV irradiation is F UV ≲ 10 3 G 0 , but only if the dust traps are located inside the photoevaporative truncation radius. The dust component of a disk will be quickly dispersed if the FUV irradiation is strong (10 4 G 0 ) or if the substructures are located outside the photoevaporation radius. We do find however, that the dust grains entrained with the photoevaporative winds may result in an absorption FUV cross section of σ ≈ 10 − 22 cm 2 at early times of evolution ( < 0.1Myr), which is enough to trigger a self-shielding effect that reduces the total mass loss rate, and slow down the disk dispersal in a negative feedback loop process.


Introduction
Protoplanetary disks are composed of the gas and dust material that orbits around newly formed stars.The evolution of isolated protoplanetary disks is often modeled by accounting for internal processes such as viscous accretion driven by magnetorotational instabilities (MRI), magneto-hydrodynamic (MHD) winds, and/or thermal photoevaporative winds due to the irradiation from the central star (Lynden-Bell & Pringle 1974;Balbus & Hawley 1991;Blandford & Payne 1982;Clarke et al. 2001;Pascucci et al. 2023).
However, growing observational evidence suggest that the irradiation from the environment also plays a significant role in the disk evolution in dense stellar regions such as the Orion Nebular Cluster, Cygnus, and Upper Sco (see Guarcello et al. 2016;Otter et al. 2021;van Terwisga & Hacar 2023, Anania et al., in prep.).In these regions, observations report particularly compact disk sizes, which points towards photoevaporative truncation, and disks that exhibit strong morphological signatures of mass loss (e.g.O'Dell et al. 1993;McCullough et al. 1995;Bally et al. 1998;Mann et al. 2014;Ballering et al. 2023).
The observational evidence is also consistent with theoretical models that predict significant mass loss rates due to the environmental radiation in a process known as external photoevaporation.(e.g.Scally & Clarke 2001;Adams et al. 2004;Anderson et al. 2013;Facchini et al. 2016;Haworth et al. 2018).
However, one of the open questions regarding the disks in these highly irradiated regions is why these disks have not yet dispersed if the mass loss rates experienced are so high, which has also been dubbed as the "proplyd lifetime problem" (Henney & O'Dell 1999).Some of the solutions proposed to this problem include taking into account different star formation events in a given region, which means that some of the disks are actually younger than the cluster itself, and also considering that stars (along with their protoplanetary disks) migrate within the cluster, and therefore they experience a varying ultra-violet (UV) irradiation during their lifetime (Winter et al. 2019a,b).Additionally, disks may be shielded from external irradiation by an optically thick envelope during the early stages of evolution, delaying the starting time of the photoevaporation process (Qiao et al. 2022;Wilhelm et al. 2023).
However, besides the survival of the gas component which is directly dispersed by external photoevaporation, it is also necessary to explain the survival of the dust component, which is observed in the millimeter continuum (Eisner et al. 2018;Otter et al. 2021).Numerical models by Sellek et al. (2020) have shown that drift is even more efficient in disks truncated by external photoevaporation, with typical depletion timescales of t depletion ≈ 2 × 10 5 yr (defined as the timescale in which 99% of the initial dust mass is lost).
In order to explain the detected millimeter fluxes and sizes, it is necessary to retain the dust for longer timescales, with dusttrapping substructures for example (Whipple 1972;Pinilla et al. 2012a).Except for very young class 0/I objects (Ohashi et al. 2023, see the recent eDisk sample), substructures such as rings and gaps appear to be a common feature in protoplanetary disks (e.g. the DSHARP sample Andrews et al. 2018), occurring even in compact disks (Zhang et al. 2023, Miley et al., in prep.).Despite this, there are currently no models that study the evolution and observability of a sub-structured disk subject to external photoevaporation.In particular, it is not clear if the solid material that is concentrated at dust traps will be able to survive the photoevaporative dispersal, or if it will be dragged along with the gas component in the thermal winds.
The goal of this paper is to characterize the evolution of the dust size distribution, the flux in the millimeter continuum, and the estimated disk size.We performed numerical simulations of disks that are subject to external photoevaporation, in cases with and without gap-(ring-)like structures.
We are interested on whether the dust traps located at the edge of gap-like structures can survive the photoevaporative dispersal process, for different UV fluxes and for different trap locations.Though we expect only the small grains to be entrained in the photoevaporative wind, the fragmentation of larger dust and uncertain growth timescale of small dust makes the outcome unclear.We also test the survival of such dust traps in the particular case of an increasing UV flux, motivated by the model Winter et al. (2019b) in which protoplanetary disks experience a variable UV irradiation during their lifetime.Finally, from our simulations we also track the distribution of the dust grains entrained with the winds, and use that information to re-estimate the effective opacities at UV wavelengths, which could result in a self-shielding effect that would protect the disk from the environmental irradiation (Facchini et al. 2016;Qiao et al. 2022).
The paper is structured as follows: In Section 2 we describe our gas and dust evolution model, which includes the growth and fragmentation of multiple grain species, the prescription for the gap-like structures, and the external photoevaporation model using the FRIED mass loss rate grid (Haworth et al. 2018;Sellek et al. 2020).In Section 3 we describe the simulations, the parameter space, and the calculation of different observable quantities.Section 4 shows our results, including a review of the difference between photoevaporating and non-photoevaporating disks, the evolution of our simulations for the different parameters, and the description of the dust content in the wind.Section 5 discusses the implications of our results and how they compare to observations of highly irradiated regions.

Disk model
We used the code DustPy 1 (Stammler & Birnstiel 2022) to simulate the gas and dust evolution of protoplanetary disks, and include the effects of external photoevaporation.

Gas evolution
The gas component evolves through viscous spreading and mass loss by external photoevaporation, through the following equation: 1 github.com/stammler/DustPy,version 0.5.6 was used for this work.where r is the radial distance to the star, Σ g is the gas surface density, and ν is the radial viscous velocity (Pringle 1981): with ν α the kinematic viscosity: which depends on the sound speed c s , the pressure scale height h g , and the turbulence parameter α (Shakura & Sunyaev 1973).
The adiabatic sound speed is defined as: where γ is the adiabatic index, T is the gas temperature, k B is the Boltzmann constant, µ is the mean molecular weight, and m p is the proton mass.
The surface density loss rate Σg,w , was derived from the FRIED grid from Haworth et al. (2018), following the implementation of Sellek et al. (2020, see their equations 5 -7).To describe it briefly, this method first predicts the expected total mass loss rate Ṁw , for a disk with surface density Σ g , around a star with mass M * , subject to an external FUV radiation field F UV , along with the photoevaporative radius r phot which is located at the interface between the optically thick and optically thin regions in the outer disk.Then, the total mass loss rate Ṁw is distributed across all the grid cells in the photoevaporative region (r ≥ r phot ), according to the material available in each grid cell.Figure 1 shows an example of the mass loss grid, with a surface density profile overlaid on top to illustrate the location of the photoevaporation radius and the regions affected by the mass loss.
This implementation assumes that the mass loss by external photoevaporation occurs only in the outer regions of the disk, and loss from the inner regions is negligible.For more details about the implementation, we refer to their original paper. 2 For simplicity, we did not include internal MHD or photoevaporative winds, which would be launched from smaller radii.

Dust dynamics
DustPy tracks the evolution of multiple grain sizes, that can grow through sticking, fragmentation into smaller species, drift towards the local pressure maximum, and diffuse according to the concentration gradient, following the model from Birnstiel et al. (2010).The corresponding advection-diffusion equation is: where Σ d is the dust surface density, d is the corresponding radial velocity, D d is the dust diffusivity, and Σd, w is the dust loss by wind entrainment.We note that we solve this equation for every dust size bin (along with the coagulation equation, see Section 2.4) and therefore all dust quantities are dependent on the grain size.
Overall, the dust dynamics can be described by its particle size a, or more specifically, by the Stokes number St, which measures the coupling between gas and dust with: with ρ s the material density of the dust, and λ mfp the mean free path of the gas molecules, where the latter is used to determine the corresponding drag regime (Epstein or Stokes I).
Given the Stokes number, the radial advection velocity of the dust is defined by: Which means that small grains (St ≪ 1) become coupled to the gas motion ν , large boulders (St ≫ 1) become decoupled from the gas, and mid-sized pebbles (St ∼ 1) drift towards the pressure maximum with a velocity of d ≈ −η K , where η = − (1/2) (h g /r) 2 dlnP/dln r measures the relative difference between the gas orbital speed and the local Keplerian speed K due to the pressure support of the gas P = Σ g c 2 s /( √ 2πh g ) (Weidenschilling 1977; Nakagawa et al. 1986;Takeuchi & Lin 2002).
We note that the effect of the pressure gradient on the dust dynamics is measured at the midplane, since large particles settle down to smaller scale heights than the gas, with (Dubrulle et al. 1995).Observational evidence of effective dust settling has recently been confirmed by observations of a handful of edge-on disks at different wavelengths (Villenave et al. 2020).The radial diffusivity is characterized by D d = ν α /(St 2 + 1), where smaller particles diffuse more easily than larger ones (Youdin & Lithwick 2007).
To include the effect of dust entrainment in the photoevaporative wind, we followed the prescription from Sellek et al. 2 The implementation of external photoevaporation in DustPy is available at: github.com/matgarate/dustpy_module_externalPhotoevaporation,and will also be included within the DustPy library package at: github.com/stammler/dustpylib. (2020), where grains smaller than the entrainment size can be lost with the wind: where Ω * = 4πh g / h 2 g + r 2 is the solid angle covered by the disk outer edge, as seen from the star.Then, the surface density loss rate for dust grains with sizes a ≤ a ent is: Σd,w (a) = ϵ(a) Σg,w , with ϵ(a) the dust-to-gas ratio of the grains in the bin size a.
In terms of disk evolution, this could lead to an outcome where the dust grains in the outer regions are either entrained at early phases during the disk lifetime, or grow and drift before they can be entrained with the wind.The study of Sellek et al. (2020) suggests that the former scenario occurs in protoplanetary disks, though differences may arise when including the full coagulation of multiple grain sizes in the model, instead of the two-population approximation from Birnstiel et al. (2012b).

Substructures and dust traps
To include the effects of substructures capable of trapping dust grains (e.g.such as the gaps formed by giant planets), we used the same approach as in Stadler et al. (2022), which implements a perturbation in the viscous α profile, that in turn creates a gaplike structure in the gas surface density profile.The perturbation in the turbulence has the shape of a Gaussian bump: where, α 0 is the turbulence base value, and A gap , r gap , and w gap are the gap amplitude, location, and width, respectively.From viscous evolution theory, a power-law gas surface density profile is scaled approximately by a factor of Σ g (r) ∝ α 0 /α(r), assuming that the disk is in steady state accretion.
The presence of a gap-like substructure then leads to the formation of a pressure maximum, where large particles (St ≳ α) can be easily trapped (Pinilla et al. 2012a, see also see Equation 7), though we note that small particles would still be able to pass through the gap by coupling with the gas component.

Grain growth
Dust growth was computed by solving the coagulation equation (Smoluchowski 1916), that accounts for the result of the collision between two grain species given their relative velocities and sizes, which can be the sticking between particles, the fragmentation of both species, and the erosion in case of a significant size difference (for more details we refer to Birnstiel et al. 2010;Stammler & Birnstiel 2022).
There are two characteristic regimes of grain growth: the drift-limited, when particles drift inward faster than they can grow, and the fragmentation-limited, when particles collide at speeds higher than the fragmentation threshold of the material frag (Brauer et al. 2008;Birnstiel et al. 2009Birnstiel et al. , 2012b)).In particular, in pressure maxima, such as the one in the outer edge of a gap, the contribution of drift to both the radial motion and the relative collision velocities between dust grains is cancelled, allowing the particles to locally accumulate and grow into larger sizes until they reach the size limit given by the fragmentation barrier.

Simulation setup
In this section we describe the initial conditions of our simulations, the parameter space explored, and the post-processing method used to derive the observable fluxes from dust grain size distribution.

Initial conditions and grid resolution
The initial surface density profile was set according to a modified version of the Lynden-Bell & Pringle (1974) self-similar solution, that include gap-like substructures that are consistent with the turbulence radial profile (see Equation 10): where M disk and r c are the initial disk mass and characteristic radius.
The gas temperature for the simulations follows a power law profile, with: where T 0 is the temperature at 1 AU.We note that this profile assumes that the heating is dominated by the passive stellar irradiation on the disk, and neglects the contribution of accretion heating, and the radiation from nearby stars.The dust distribution is initialized assuming a uniform constant dust-to-gas ratio, with Σ d = ϵ 0 Σ g , following the ISM size distribution from (Mathis et al. 1977), which goes from 0.5 µm to an initial maximum grain size of a 0 = 1 µm.
The radial grid was set to be linearly spaced with r 1/2 , and going from 2.5 AU to 500 AU, with n r = 200 grid cells.The FRIED grid is tabulated out to 400 AU, however, in our models the photoevaporative truncation radius is always smaller than the outer radial grid boundary of the FRIED grid and no extrapolation in the 400-500 AU range is needed.The grid for the dust size distribution was set with a logarithmic spacing, going from approx.0.5 µm to 50 cm with n m = 127 grid cells, such there are 7 grid cells for each order of magnitude in mass in order to reliably solve the grain coagulation (Ohtsuki et al. 1990;Dr ążkowska et al. 2014;Stammler & Birnstiel 2022).We evolved the simulations from time t = 0 up to 5 Myr, though some disks may fully disperse earlier due to the influence of external photoevaporation.

Parameter space
In this work, our main focus is to explore on how the gap presence and its location affects the retention of solids in protoplanetary disks that are subject to different FUV radiation environments.For the fiducial model we took a 1 M ⊙ mass star, surrounded by a disk with mass of M disk = 0.1M * , and an initial characteristic radius of r c = 90 AU.We started our study with an initial comparison against disks without photoevaporation, to have an overview of the key aspects of the disk evolution.
For the main parameter space exploration, gap-like substructures can be located at 1/3 r c (inner trap) or 2/3 r c (outer trap), or completely absent.(A gap = 0, no traps).With r c = 90 AU, this means that the gaps are located either at 30 AU or 60 AU in our simulations.Bae et al. (2022) compiled data of disks with observed substructures, showing that most of the rings have been found between 20-60 AU (their histogram in Fig. 3d), as assumed in this work.Some disks have substructures up to 100-200 AU, which are outside the photoevaporation radius in our models.The disk can be subject to external photoevaporation due to far ultra-violet (FUV) fluxes of 10 2 G 03 (weak), 10 3 G 0 (medium), or 10 4 G 0 (strong) (where G 0 corresponds to the local interstellar radiation field Habing 1968).We note that in low mass star-forming regions such as Taurus/Lupus the external FUV radiation field is of order 1-100 G 0 .In more massive stellar clusters such as the Orion Nebular Cluster the FUV radiation field strength that disks are exposed to ranges from 1 − 10 7 G 0 .Previous work has suggested that the most common experienced UV environment is around 10 3 G 0 (Fatuzzo & Adams 2008;Winter et al. 2020).Additionally, we also study how the disk evolution changes for a star with lower mass (correspondingly with lower disk temperature), a disk that is initially more compact, and a disk that is subject to a FUV radiation field that increases with time (in order to mimic, for example the effect of migration within a stellar cluster, or the clearing of the original molecular cloud Winter et al. 2019b;Qiao et al. 2022).The complete parameter space and the additional disk physical parameters can be found in Table 1.Grain material density 1.67 g cm −3

Fluxes in the millimeter continuum
To compare our models with observations, we obtain the grain size distribution from the dust evolution simulations with DustPy, which is then used to compute the intensity profile and the total flux in the millimeter continuum, at λ = 1.3 mm.We assume that the dust grains follow the opacity model from Ricci et al. (2010), and are composed of water ice, carbon and silicates (Zubko et al. 1996;Draine 2003;Warren & Brandt 2008).
With the opacity profile k ν (a), with ν the frequency, we can calculate the optical depth at every radius as: the intensity profile with: and the total flux as F ν = I ν dΩ, where dΩ is the differential of the solid angle covered by the disk in the sky, and B ν (T ) is the emission of a black body with temperature T .For the purpose of this work, we assumed our disks to be at a distance of 400 pc, which is the approximate distance to the Orion Nebular Cluster (ONC, Hirota et al. 2007;Kounkel et al. 2017), and to account for observational limitations, we also assume a beam size of 40 mas and a sensitivity threshold of 0.1 mJy beam −1 (this is used to calculate the measured disk size r 90% , and the intensity profile I ν ).While for the gas and dust evolution, both of them have the same temperature (Eq.12), for the calculations of the millimetre fluxes, we assume that there is a background temperature of T b = 20 K, to account for the irradiation from massive stars in the interstellar space, thus:

Results
In this section we show how the mass evolution, dust size distribution, and observable properties such as the intensity profile, total flux and the observable dust disk size change due to the effect of external photoevaporation and substructures within our parameter space.

Fiducial comparison
We begin our comparison between simulations with and without external photoevaporation, in order to understand the main features of each evolution pathway.The effect of photoevaporation on the gas component is straightforward (Figure 2, top panel).The gas mass decreases faster for disks undergoing external photoevaporation than in the viscous evolution counterparts, and (for our model) the presence of substructures have no significant impact on the overall gas evolution.Meanwhile, the dust component is distinctively affected by both the effects of photoevaporation and dust traps.Figure 3 shows the morphology of the dust size distribution at relevant times, where we see how photoevaporation truncates the outer disk, and how dust traps retain higher concentrations of large grains.Disks subject to external photoevaporation have lower dust masses than both of their viscous evolution counterparts (Figure 2, middle panel), and the mass loss in (externally) photoevaporating disks can be described in two stages: a wind dominated loss, and a drift dominated loss (Figure 4).During the first stage of disk evolution, the dust grains in the outer regions are still small and can be easily entrained with the photoevaporative winds.This leads to a decrease in the dust mass from the initial 300 M ⊕ to 80 M ⊕ during the first 0.1 Myrs.In contrast, the dust mass in viscous evolution simulations only decreases to 250 M ⊕ during the same period of time.
From dust evolution theory, the lifetime of the disk (in terms of the dust component) is on the same order of magnitude as the drift and growth timescales of dust particles at the disk outer edge (Birnstiel et al. 2012a;Powell et al. 2017).Because photoevaporation removes all the material from the disk outer regions, effectively truncating the disk size, the remaining dust component will drift faster towards the star, in comparison with the viscous counterparts.
Figure 4 shows that in the simulation with photoevaporation and without dust traps, the dust loss rate is first dominated by M dust /M gas Fig. 2. Gas and dust mass evolution (top and middle panels), and the global dust-to-gas ratio (bottom panel).The plot shows the evolution of a disk with and without the influence of FUV external photoevaporation (black vs. gray lines), for disks with and without a gap-like substructure (solid vs. dotted lines).For this comparison, the disk is around a M ⊙ mass star with an initial size of r c = 90 AU.The gap substructures are located at r gap = 1/3r c , and the irradiation field is wind entrainment (t ≲ 0.1 Myr), and quickly becomes dominated by drift (t ≳ 0.1 Myr).For this simulation, the disk has lost 99% of the initial dust mass by t depletion ≈ 0.4 Myrs.Up until this point, our results agree with the work of Sellek et al. (2020, see their Fig.5), and we refer to their paper for a detailed analysis of the evolution of disks without substructures.
It is during the drift dominated stage that the presence of a dust trap becomes relevant, greatly delaying the depletion of the dust component.In this simulation with external photoevaporation and dust traps, the total mass lost by drift becomes comparable to that lost by winds after t = 1 Myr, and the disk loses the 99% of the initial dust mass only by t depletion = 2.3 Myrs.This depletion timescale exceeds the values reported in Sellek et al. (2020) by an order of magnitude, indicating that the presence of dust traps may explain why disks in dense star clusters subject ] Fig. 3. Dust size distribution at 0.1, 1, and 3 Myr (from left to right), for the simulations with/without photoevaporation, and with/without dust traps.The solid lines indicate the estimated fragmentation and drift growth limits (magenta and cyan, respectively).For reference, we show the grain sizes that correspond to St = 0.1 with a dashed white line, and the photoevaporation radius with a dashed red line.Note that we plot σ dust , which has surface density units, but it is normalized by the logarithmic bin size (see Birnstiel et al. 2010).to high F UV fluxes are still detectable in observations (Guarcello et al. 2016;Eisner et al. 2018;Otter et al. 2021).
In terms of the global dust-to-gas mass ratio (Figure 2, bottom panel) it is interesting to note that this is always below the initial ϵ 0 = 0.01, despite the gas mass loss in disks with external photo-evaporation.In particular, the simulation without dust traps under photoevaporation displays a remarkably lower dustto-gas ratio than its viscous counterpart, which is due to the re-duced drift timescale by the photoevaporative truncation of the disk.
The effect of photoevaporation and the dust traps can also be seen in the derived observable quantities, i.e. the intensity profiles, the flux, and the dust disk radius that is usually assumed as the radius enclosing 90% of the total millimeter flux (obtained following Equation 14 corresponding fluxes, and the radii containing 90% of the dust emission at specific times in Table 2. From the intensity profiles (top panel, shown at 1 Myr) we see how the disk subject to external photoevaporation are truncated at r ≈ 60 AU − 70 AU, in contrast with their viscous evolution counterparts, which display a more extended emission.As expected from the difference in dust masses, the disks with the gap-substructures are brighter than their smooth counterparts, featuring a ring-like emission at r ≈ 35 AU.We note that, since the dust trap is located well inside the truncation radius, the morphology of the bright ring is not affected by external photoevaporation.
The evolution of the total flux F ν (middle panel) follows the trend of the total dust mass.Here, the simulations with photoevaporative loss show fluxes that range between 0.01-15 mJy after t = 1 Myrs, and in particular the simulation with a dust trap manages to display an emission higher than 1 mJy until 4 Myrs.The value of r 90% (bottom panel) reflects the extension of the intensity profile until the sharp drop below the sensitivity limit (0.1 mJy beam −1 for this plot).For the disks subject to external photoevaporation the disk size progressively shrinks, first until the truncation radius within the first 0.1 Myr due to the initial dust removal by winds.Then, if no dust traps are present, the disk continues to shrink until it disappears at t ≈ 1 Myr (i.e., the emission is fully below the sensitivity limit) or, if substructures are present, the disk only shrinks until the location of the dust ring.
We also note that the viscous evolution simulations first experience an expansion until reaching r 90% = 150 AU at t = 0.5 Myr, and then contract as the dust drifts inwards, as the emission from outer regions of the disk falls below the sensitivity limit.At 3 Myr the dust disk sizes in the viscous simulations either reach the dust trap location or drop sharply (for the case without substructures).

Parameter Space: UV field and Trap location
In this section, we study how the presence of dust traps affects the dust mass and observable properties of the disk depending on the external FUV flux (Figure 6 and Figure 7), and also list the dust depletion timescale (i.e. the age of the disk when 99% of the initial dust mass has been lost) covering the parameter space of Table 3.The table also indicates whether a dust trap  3 mm continuum at a distance of 400 pc.Bottom: Time evolution of the r 90% enclosing 90% (black) of the disk emission in the continuum above the sensitivity limit (0.1 mJy beam −1 ), and the photoevaporative radius derived from the FRIED grid mass loss profile (red dotted lines, see Sellek et al. 2020).We note that towards the end of the simulation the calculation of the photoevaporative radius suffers from numerical effects due to the discrete nature of the FRIED grid.
was dispersed by the photoevaporative winds within the 5 Myr of simulation time.
For the strongest FUV flux (10 4 G 0 ), the mass loss due to photoevaporation is so intense, that it renders existing substructures irrelevant in terms of observability.From the evolution of r 90% and the photoevaporative radius we can see how photoevaporation progressively shrinks the disk, until it completely disappears by 1 Myr, approximately at the same time that the emission in the millimeter continuum drops.We see that if the substructure is located in the inner regions (r gap = 30 AU, or 1/3 r c in our simulation), there is only a minor delay of approx 0.2 Myr until the flux in the millimeter continuum drops.
For the medium FUV flux (10 3 G 0 ), we already saw in the previous section that an inner dust trap can extend the disk lifetime by 2 Myr, and that the disk r 90% size converges to the location of the dust trap, as it also occurs in models with dust traps However, we also find that in the case where the substructure is located further out (r gap = 60 AU, or 2/3 r c ), the disk would still be dispersed by photoevaporation, once the photoevaporative radius reaches the location of the dust trap.This can be seen in Figure 7 (mid-bottom panel, dashed line), where the disk size initially matches the location of the outer dust trap, and quickly disperses once the disk r 90% becomes comparable with the photoevaporative radius.In terms of the dispersal timescale and continuum fluxes, there is almost no difference between having an outer dust trap or no dust trap at all.In the case with a weak FUV field (10 2 G 0 ), both the inner and outer dust trap remain inside the photoevaporation radius (as seen from the value of r 90% ), which means that the size and flux of the disk will be dominated by the material trapped at local the pressure maximum.In comparison with the medium and strong photoevaporative regimes, disks in the weak photoevaporative environment have longer lifetimes of t depletion = 3 Myr and t depletion ≳ 5 Myr, for the cases with inner and outer substructures respectively.For this regime, the lifetime of the disk with the outer substructure is longer than the one with the inner substructure, which we infer to be due to the longer drift timescales from the outer dust trap to the star, and the longer diffusion timescales across the gap in the outer regions.
From these three photoevaporative regimes we conclude that the lifetime of the disk is dominated by the outermost dust trap, provided that this is located well inside the photoevaporative radius, and that the dust trap located in the photoevaporative regions (even marginally) would have little to no effect on the disk survival timescale.It might seem surprising that disks with fully formed dust traps are dispersed, since photoevaporative winds can only carry away small grains.However, because fragmentation of large grains continuouslyreplenishes the population of micron-sized particles, there is always a non-negligible fraction of solids being carried away.We will revisit the dust size distribution of entrained grain in Section 4.5.

Low mass stars and compact disks
We perform two additional simulations to study whether a dust trap would survive in a disk around a low mass star (M * = 0.3 M ⊙ ), and if an initially more compact disk size (r c = 60 AU) would affect our observational predictions.For this comparison, we use a radiation field of F UV = 10 3 G 0 , and fix the dust trap location at 30 AU.
In Figure 8 we see that the disk around a 0.3 M ⊙ mass star is dispersed faster than the one around a sun-like star, and that even the presence of the dust trap cannot prevent the sharp decrease in the millimeter continuum flux.This faster dispersal oc- curs because the gravitational potential is weaker around a low mass star, which means that the gas and solid can be removed more easily by the UV irradiation from the environment.Consequently, the photoevaporation radius can reach further into the inner regions of the disk, all the way down to r phot ≈ 20 AU, and disperse the dust trap that was located at 35 AU.In order for the dust component to survive in a disk around low mass star, the trap should be located further inward (r gap ≲ 10 AU), provided that the UV field is on the order of F UV = 10 3 G 0 .Disks in regions with lower irradiation could still display dust traps at larger radii.
For the case of a disk that is initially more compact, we do not see any significant differences in the disk evolution in terms of the observed size r 90% or the flux in the millimeter continuum (see Figure 9).From this plot, we can expect for disk sizes in photoevaporative regions to be determined by their outermost dust trap, independently of the initial extension.et al. (2019b) showed some of the disks in high irradiation regions, could have actually formed and migrated from lower irradiation environments.This would explain why these objects are still observable despite the short dispersal timescales associated with the high irradiation.

Winter
We conduct a simple simulation with our fiducial parameters (M * = 1 M * , r c = 90 AU, r gap = 30 AU) in which we gradually increase the F UV irradiation with the following function:   16).As in Figure 3, the fragmentation limit is marked in magenta, the drift limit in cyan, and the photoevaporation radius as a vertical dashed-red line.
where F UV,0 = 10 G 0 and F UV,max = 10 4 G 0 are the lower and upper limits of the FRIED grid (Haworth et al. 2018).The exponent of 3/2 is meant to represent an UV flux that increases rapidly as the disk approaches the high irradiation regions, though the exact shape would depend on variables such as the disk trajectory, the distribution of bright massive stars within the cluster (Winter et al. 2019b;Qiao et al. 2022;Wilhelm et al. 2023), and the variation in their luminosity with time (Kunitomo et al. 2021).
In Figure 10 we show the evolution UV flux and the mass loss rate over time up to 1 Myr.16).
From the dust distribution shown in Figure 11 we see how the radial extension of the dust component gets progressively smaller as time passes, and in particular, we see how the dust trap completely vanishes by 1 Myr once the UV flux reaches its peak.We note that the dust trap dispersal occurs despite the fact that dust grains have already grown larger into millimeter-tocentimeter sized pebbles (Figure 11, middle panel).While these particles are not easily entrained by the wind, they are still being indirectly depleted, since they continuously fragment into the smaller size grains that are directly removed by photoevaporation.
This phenomenon can also be seen in the spike in dust loss rate at 0.7 Myr (Figure 10, bottom panel), which coincides whith the photoevaporative radius when it catches up with the location of the dust trap, and in the millimeter continuum (Figure 12), when the flux sharply drops.From our results in Table 3 for the high UV fluxes (10 4 G 0 ) we can expect the remaining dust component to disperse in timescales of 0.2 Myr or less.This implies that even if dust traps manage to survive during an early low irradiation, for example if the disk was initially shielded from UV irradiation (Qiao et al. 2022;Wilhelm et al. 2023), farther away from the dense regions of the cluster (Winter et al. 2019b), or surrounded by intermediate mass stars (M ≲ 3M ⊙ ) in early evolutionary stages (when the UV emission is lower Kunitomo et al. 2021), once the FUV flux increases the dust component in the dust traps will be quickly dispersed along with the gas.Size ( m) ] 10 1 10 2 r (AU) 10 0 10 1 10 2 Size ( m) ] Fig. 13.Grain size distribution (integrated in time) of the dust grains lost with the photoevaporative flow for the disks subject to low, medium and high F UV environments, at t = 1 Myr.The disk's parameters are a star of 1M ⊙ , initial size of r c = 90 AU, and a gap substructure at 30 AU.The location of the gap r gap is indicated with a dotted white line.

Dust distribution of the lost grains
Photoevaporative winds can remove the small grains entrained with the gas flow from the disk outer regions (defined through the photoevaporative radius).While tracking the subsequent spatial evolution of these grains goes beyond the scope of this paper, we can still record of the mass and size distribution of lost dust component, and compare between the different irradiation regimes.
Figure 13 shows the time integrated distribution of all the dust grains that have been lost by entrainment with the wind up to 1 Myr, as a function of the grain size and the launching location, that is: where we see how the maximum grain size entrained in the winds and the launching region depend on the irradiation from the environment.In particular, we see how for the weakest irradiation field (10 2 G 0 ) only grains smaller than 10 µm can be removed from regions beyond 100 AU, while for the strongest irradiation field (10 4 G 0 ) even grains of sizes up to a ≈ 100 µm can be dragged along by the winds, with a noticeable increase in the dust removal at r ≈ 30−40 AU where the dust trap is located.
To obtain a broader overview of the lost grain distribution, we show in Figure 14 the cumulative size distribution for the parameter space shown in Section 4.2, i.e.: which shows that approximately 70% − 80% of the lost mass of solids comes from the sub-micron size particles.Overall, we find that the location of the dust trap has little to no impact in the size distribution of the lost grains (specially in the sub-micron size range).For the weak irradiation environment (10 2 G 0 ) the dust traps are located well inside the photoevaporation radius and do not contribute to the dust content in the wind, and for the strong photoevaporative environment (10 4 G 0 ) there is an increase of 15 M ⊕ in the micron range size range for configuration with the inner dust trap respect to the case without dust trap.The dust distribution within the photoevaporative winds should also, in principle, determine the effective cross section and opacity of the disk to the environment FUV irradiation (Facchini et al. 2016).If the dust content and effective opacity are high, then the disk could be self-shielded (Qiao et al. 2022;Wilhelm et al. 2023)  we calculate what would be the self-consistent cross section for FUV photons per gas molecule with: where ϵ w (a) is the dust-to-gas ratio of the dust species with size a in the wind, κ UV (a) is the absorption opacity at λ = 0.1 µm, µ = 2.3 is the mean molecular weight, and m p is the proton mass (see also Facchini et al. 2016, Eq. 23).In ?? we show the dust-togas ratio in the photoevaporative wind, and the respective FUV cross section per gas molecule (using the Ricci opacities as in the results in Sect.4).We find that the total dust-to-gas ratio remains almost constant in the first ∼ 0.05-0.1 Myr of evolution and close to the initial value of 10 −2 for all the irradiation regimes and trap configurations.Once the dust has had a chance to grow, the dust-togas ratio sharply decreases with time as the dust loss rates due to photoevaporation decreases over time (see Fig. 15) Following the same trend, the effective cross section for UV wavelengths is on the order of σ FUV ≈ 2 − 2 × 10 −22 cm 2 .This is calculated assuming We also perform the calculation of the FUV cross section assuming an ISM grain size distribution (Mathis et al. 1977) for comparison (with σ ISM ≈ 2.5 × 10 −22 cm 2 ), and find that all our simulations fall below this value, specially after grains have growth after ∼0.02-0.1 Myr.Though grain growth is expected to reduce the effective absorption cross section at FUV wavelengths (Facchini et al. 2016) since larger grains have lower absorption opacities, the decrease of σ FUV overtime seems to be mostly dominated by the decrease of the dust-to-gas ratio in the wind.
In section Section 5.2 we compare the values found for the FUV cross section with previous studies, and discuss if selfshielding by the grains entrained in the wind could be important for the global disk evolution and dispersal process.

Discussion
5.1.Can dust traps explain the observations of highly irradiated regions?
Observations in the millimeter continuum by Eisner et al. (2018); Otter et al. (2021) of the ONC and OMC1 star-forming regions show protoplanetary disks with fluxes on the order of 0.1 mJy to 10 mJy at 0.85, 1.3, and 3 mm wavelengths, as well as a lack of disks sizes larger than 75 AU (Otter et al. 2021) which could be explained by photoevaporative truncation due to external photoevaporation.
However, the work of Sellek et al. (2020) showed that the lifetime of the dust component is significantly shorter in disks undergoing external photoevaporation, than in disks undergoing regular viscous evolution (see Section 4.1), with depletion timescales that are on the order of 0.1 Myrs.This means that the survival of the dust component in observations could not be explained even in medium radiation environments without an additional process that prevents the dust loss.
Since the study of Sellek et al. (2020) does not consider the presence of substructures in protoplanetary disks, which currently are known to be a common feature (e.g.DSHARP sample Andrews et al. 2018), we test whether the presence of dust traps could be help to explain the lifetime of the dust component in photoevaporative disks.Though it might seem evident that substructures should contribute to increase the disk lifetime(e.g.Pinilla et al. 2012a,b), it is not clear whether the dust grains will be able to resist the drag force from the gas in the photoevaporative winds, specially considering that fragmentation continuously replenishes the population of the small grains that are easily entrained with gas.
We find in this work that a gap-like substructure can significantly increase the dust component lifetime in weak radiation environments (10 2 G 0 ) to more than 5 Myr, and that in medium radiation environments (10 3 G 0 ) it can increase the disk lifetime to approx. 2 Myrs, but only if the dust trap is located well inside the photoevaporation radius, which in our parameter space was approximately between 50 and 75 AU.The simulations where the dust traps were located outside the photoevaporation radius were dispersed, as the small grains were replenished by fragmentation at a faster rate than the spatial evolution of the photoevaporation front, resulting in the dispersal of the dust traps.
For extreme photoevaporation environments (10 4 G 0 ) we observe that substructures have no significant effect on the survival timescale of the protoplanetary disk solid component.The photoevaporation front quickly truncates the whole disk from the outside-in, dragging all the solid grains along with the wind.Because this strong irradiation regime is precisely the one that affects the disks in dense stellar regions such as the core of the ONC, we infer that dust traps alone would not be able to explain the millimeter emission in this type of environment.
We also consider the possibility that photoevaporation might be delayed, with the environmental FUV flux increasing over time instead of acting at full strength from the very beginning of the simulation.This scenario would be consistent with migration across the stellar cluster (Winter et al. 2019b), or shielding of the disk by the primordial envelope (Qiao et al. 2022), however in this case we still find that the material in the dust traps is quickly dispersed once the photoevaporation front reaches the location of the pressure maximum (see Section 4.4).
We conclude that the presence of dust traps can increase the disk lifetime in weak and mild photoevaporation environments, which in the case of disks without substructures is limited by the drift timescale at the truncation radius, but cannot help to retain the dust component beyond the lifetime of the gas component, as the grains are efficiently dragged along with the photoevaporative winds.Therefore, while dust traps seem necessary to prevent the dust from draining too quickly by drift, these are not sufficient to explain the survival of the disks in the millimeter continuum, and models that extend simultaneously the lifetime of the gas and dust component such as those accounting for shielding, migration, or multiple star formation events (Qiao et al. 2022;Wilhelm et al. 2023;Winter et al. 2019a,b) are better suited to explain the effective actual disk lifetime.Other mechanisms that could contribute to the disk long-term survival are the luminosity evolution of bright B stars, which require approximately 1 Myr from their formation to reach their peak brightness at UV wavelengths (see Kunitomo et al. 2021, Fig. 3), and dominate the irradiation field in regions such as Upper Sco (Anania et al., in prep.),It is necessary to conduct high resolution observations in the millimeter continuum of disks that can be subject to influenced by high environmental irradiation, and identify first if substructures are present, and second if their presence (or absence) would be consistent with the age of the dust component.One promising target for this study would be the disk ISO-Oph2 (Casassus et al. 2023), which is in the proximity of the B star HD147889, and shows signatures of being heated by its bright neighbour.

Self-shielding by the grains entrained in the wind
The calculations performed in this paper use the FRIED grid from Haworth et al. (2018), which assumes that photoevaporative winds are depleted in dust, with a low dust-to-gas ratio of ϵ w = 3×10 −4 and a FUV cross section of σ FUV = 2.7×10 −23 cm 2 .In Fig. 15 we recalculated the FUV cross section using the dustto-gas ratio measured directly from the simulation lost material, and found them to be higher (σ FUV ≈ 10 −22 cm −2 ) at early times of evolution (∼0.05 Myr) than the values used by Haworth et al. (2018).This means that, in order to be self-consistent with the dust content in the wind, the effective mass loss rates should be lower than those used for this paper (Haworth et al. 2023) This self-shielding effect is a negative feedback loop which would moderate the strength of the photoevaporative dispersal.As such we could expect mass loss to proceed at a slower rate, and prolong the disk lifetime in both the gas and dust component, though a dedicated study in which the mass loss and cross section are computed simultaneously in run time would be necessary to confirm this theory.
Another important effect to account when accounting for self-shielding would be the spatial evolution of the dust and gas after these are launched from the disk surface (Paine et al., in prep.).If the ejected material disperses quickly (as suggested in our results in Fig. 15), it will not be able to shield the disk from the external irradiation sources since the cross section sharply drops after the first few ∼0.02-0.1 Myr of disk evolution.Therefore, in order to determine if self-shielding can extend the disk lifetime or not, it is also necessary to study whether the dust and gas material lifted by the photoevaporative winds can remain around their parent protoplanetary disk for longer timescales, on the scale of 1 Myr or more.
We note that in this work we assumed the opacity values of Ricci et al. (2010) to calculate the FUV cross sections, while previous calculations of Facchini et al. (2016, for example), used the optical constants from Li & Greenberg (1997) which leads to the difference in the reported for the ISM cross section, where our estimation is lower approximately by a factor of 4. We infer that our cross sections would increase by a factor of a few if we had used the same optical constants and grain composition.
Previous studies have shown that shielding can reduce the effect of external photoevaporation on the disk (e.g.Qiao et al. 2022), however its effect is usually constrained to the early stages of disk evolution.In contrast, if the self-shielding by the wind-entrained grains is proven to be effective, it could actually contribute to reduce the mass loss in high irradiation environments for extended periods of time, specially in disks with multiple substructures.
We highlight the importance of developing self-consistent models that couple dust trapping, photoevaporation driven by external irradiation, and the dust content in the ejected winds, since these ingredients do interact with each other and change the course of the disk evolution, as shown, for example, in Owen & Lin (2023) where it was proposed that circumstellar disks in extreme environments such as the galactic centre could survive photoevaporative dispersal.

Summary
In this work, we studied the evolution of the dust component and emission in the millimeter continuum of protoplanetary disks subject to high environmental FUV irradiation, accounting for the presence of gap-like substructures that act as dust traps.As in Sellek et al. (2020), we also find that dust drift is more efficient in disks subject to external photoevaporation than in standard viscous disks, where the lifetime of the dust component can be as short as a few 0.1 Myr if dust traps are not present.
In weak irradiation environments of F UV = 10 2 G 0 , the presence of dust traps does prevent drift of dust grains into the star and allow the disk to survive for several Myrs.In irradiation environments of 10 3 G 0 , dust traps need to be inside the photoevaporative truncation radius to extend the disk lifetime (from 0.3 to 2.3 Myr for our choice of parameters, see Section 4.2), while dust traps located outside the truncation radius are dispersed with the photoevaporative winds and do not extend the lifetime of the dust component significantly.Finally, in extreme irradiation environments of 10 4 G 0 all the dust traps are dispersed as the photoevaporation front clears the entire disk from the outside-in.
Though dust traps are a necessary ingredient to explain the survival of the dust and the observed millimeter emissions in highly irradiated environments, specially considering that drift is specially efficient in disks truncated by external photoevaporation (Sellek et al. 2020), they are not enough to explain why the objects found at the core of dense stellar regions such as the ONC or Cyg OB2 (Guarcello et al. 2016;Eisner et al. 2018;Otter et al. 2021) have not yet dispersed.
Overall, it seems more likely that the disks observed in these regions have only recently began to experience the high irradiation and extreme mass loss rates measured in observations.These objects could have been subject to a much lower irradiation earlier in their evolution due to shielding, migration from the regions with lower stellar densities, or they could have been born in a more recent star formation event than the rest of the cluster (Winter et al. 2019a,b;Qiao et al. 2022).
On that line, we do find that the dust content entrained with the photoevaporative winds might be enough to increase the cross section for FUV photons to σ FUV ≈ 10 −22 cm 2 at 0.1 µm wavelengths at early times of evolution, partially shielding the protoplanetary disk from external irradiation, and decreasing the total mass loss rate.However, it is yet unclear whether the dusty material lifted with the winds will stay in the proximity of the parent disk, effectively acting as a shield, or if it will quickly disperse leaving the disk surface bare to the environmental irradiation, and further studies are necessary.
In conclusion, dust traps are necessary but not sufficient to explain the survival of the dust component and continuum emission in high photoevaporative regions, and self-shielding by the dust entrained with the wind may contribute to reduce the extreme mass loss rates, but only if the dusty material can remain close enough to the parent disk to absorb the irradiation from the environment.

Fig. 1 .
Fig.1.Example mass loss rate grid for a disk orbiting a 1 M ⊙ star subject to an irradiation of F UV = 10 3 G 0 .The grid determines the gas mass loss Ṁw as a function of the gas surface density Σ g and radii r, using the FRIED grid fromHaworth et al. (2018) and the implementation ofSellek et al. (2020).The solid white line, shows an example gas surface density profile, the dashed vertical line indicates the photoevaporative truncation radius r phot , and the solid red line indicates the regions which are subject to external photoevaporation.The photoevaporation radius evolves with the simulation, and corresponds to the maximum of the mass loss rate along the current surface density profile.

Fig. 4 .
Fig.4.Top: Evolution of the dust loss rate over time for a photoevaporating disk (with F UV = 10 3 G 0 ), with and without dust traps (solid vs. dotted lines).The figure shows the distinction between the dust loss rate due to drift into the star (red), and due to entrainment with the photoevaporative wind (blue).Bottom: Cumulative fraction of dust lost (relative to the initial dust mass) due to winds and drift.

Fig. 5 .
Fig. 5. Top: Intensity profile for the fiducial simulations at λ = 1.3 mm, assuming a distance of d = 400 pc, and a beam size of 40 mas (16 AU, shown as a horizontal blue line), for a time of t = 1 Myr.Middle: Evolution of the disk flux at λ = 1.3 mm.Bottom: Evolution of the radius enclosing 90% of the disk continuum emission, assuming a sensitivity threshold of 0.1 mJy beam −1 .

Fig. 6 .
Fig.6.Top: Mass evolution for our parameter space in UV fluxes and trap location (or absence).Bottom: Fraction of dust mass lost.The lines distinguish between the total dust lost (black), the fraction lost by drift at into the star (blue), and lost by entrainment with the photoevaporative wind (red).

Fig. 7 .
Fig. 7. Evolution of the disk observables, for different FUV fields and trap locations (1/3r c , 2/3r c or none).Top: Evolution disk flux in the 1.3 mm continuum at a distance of 400 pc.Bottom: Time evolution of the r 90% enclosing 90% (black) of the disk emission in the continuum above the sensitivity limit (0.1 mJy beam −1 ), and the photoevaporative radius derived from the FRIED grid mass loss profile (red dotted lines, seeSellek et al. 2020).We note that towards the end of the simulation the calculation of the photoevaporative radius suffers from numerical effects due to the discrete nature of the FRIED grid.

Fig. 8 .
Fig. 8. Evolution of the disk continuum flux at λ = 1.3 mm and the r 90% (assuming a sensitivity threshold of 0.1 mJy beam −1 ) for two different stellar masses, and including the gap at r c = 30 AU.The photoevaporative radius is shown for comparison.The UV flux is 10 3 G 0 .

Fig. 11 .
Fig. 11.Dust size distribution at 0.1, 0.5, and 1 Myr for the simulation with an increasing UV Flux (Equation16).As in Figure3, the fragmentation limit is marked in magenta, the drift limit in cyan, and the photoevaporation radius as a vertical dashed-red line.

Fig. 12 .
Fig.12.Same as Figure8, comparing the simulation with constant UV Flux of 10 3 G 0 and the one with increasing UV flux (Equation16).

Fig. 14 .
Fig.14.Cumulative size distribution of the dust grains lost up to a time of 1 Myr for the simulations shown in Section 4.2, for the different UV fluxes, and trap properties.The line styles represent the disks with an inner trap, an outer trap, and no traps at all (solid, dashed, and dotted lines respectively).For the disks with low irradiation environments (10 2 G 0 ) the dust traps have no effect on the lost dust distribution and the three lines overlap.

Fig. 15 .
Fig.15.Top: dust-to-gas ratio of the material removed by the photoevaporative winds for different irradiation environments and dust trap configurations.Bottom: FUV cross section per gas molecule of the material removed by stellar winds, considering the dust-to-gas ratios from the top panel.

Table 1 .
Parameter space.Fiducial values are highlighted in boldface.

Table 3 .
Dust depletion timescales."Dispersed" refers to dust traps that have vanished in the simulation due to external photoevaporation.
, leading to a negative feedback loop where photoevaporation regulates itself.As a first step to study this process