© The Authors 2024

1 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 magneto-rotational 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 pho-toevaporation. (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⁵ 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 dust-trapping 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 Sect. 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 Sect. 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.

2 Disk model

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

2.1 Gas evolution

The gas component evolves through viscous spreading and mass loss by external photoevaporation, through the following equation: $$\frac{\partial }{\partial t}{\text{Σ}}_{\text{g}}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{\text{Σ}}_{\text{g}}{v}_v\right)=-{\dot{\text{Σ}}}_{\text{g},\text{w}},$$(1)

where r is the radial distance to the star, Σ_g is the gas surface density, and υ_v is the radial viscous velocity (Pringle 1981): $$v_v=-\frac{3}{{\text{Σ}}_{\text{g}}\sqrt{r}}\frac{\partial }{\partial r}\left(v_{\alpha }{\text{Σ}}_{\text{g}}\sqrt{r}\right),$$(2)

with v_α the kinematic viscosity: $$v_{\alpha }=\alpha c_{\text{s}}{h}_{\text{g}}\text{, }$$(3)

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: $$c_{\text{s}}=\sqrt{\gamma \frac{k_{\text{B}}T}{\mu m_{\text{p}}}},$$(4)

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 ${\dot{\Sigma }}_{\text{g},\text{w}}$, was derived from the FRIED grid from Haworth et al. (2018), following the implementation of Sellek et al. (2020, see their Eqs. (5)–(7)). To describe it briefly, this method first predicts the expected total mass loss rate ${\dot{M}}_{\text{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 ${\dot{M}}_{\text{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². For simplicity, we did not include internal MHD or photoevaporative winds, which would be launched from smaller radii.

Fig. 1 Example mass loss rate grid for a disk orbiting a 1 M⊙ star subject to an irradiation of FUV = 103 G0. The grid determines the gas mass loss ${\dot{M}}_{\text{w}}$ as a function of the gas surface density Σg and radii r, using the FRIED grid from Haworth et al. (2018) and the implementation of Sellek et al. (2020). The solid white line, shows an example gas surface density profile, the dashed vertical line indicates the photoevaporative truncation radius rphot, 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.

2.2 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: $$\frac{\partial }{\partial t}{\Sigma }_{\text{d}}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{\Sigma }_{\text{d}}{v}_{\text{d}}\right)-\frac{\partial }{\partial r}\left(r{D}_{\text{d}}{\Sigma }_{\text{g}}\frac{\partial }{\partial r}\left(\frac{{\Sigma }_{\text{d}}}{{\Sigma }_{\text{g}}}\right)\right)=-{\dot{\Sigma }}_{\text{d},\text{w}},$$(5)

where Σ_d is the dust surface density, υ_d is the corresponding radial velocity, D_d is the dust diffusivity, and ${\dot{\Sigma }}_{\text{d},\text{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 Sect. 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: $$\text{St}=\frac{\pi }{2}\frac{a{\rho }_{\text{s}}}{{\Sigma }_g}\cdot \{\begin{array}{ll}1\hfill & {\lambda }_{\text{mfp}}/a\ge 4/9-\text{ Epstein},\hfill \\ \frac{4}{9}\frac{a}{{\lambda }_{\text{mfp}}}\hfill & {\lambda }_{\text{mfp}}/a<4/9-\text{ Stokes }I,\hfill \end{array}$$(6)

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: $$v_{\text{d}}=\frac{1}{1+{\text{St}}^2}{v}_v-\frac{2\text{St}}{1+{\text{St}}^2}\eta v_k.$$(7)

Which means that small grains (St ≪ 1) become coupled to the gas motion υ_v, 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 η = − (l/2> (h_g/r)² din P/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={\text{Σ}}_{\text{g}}{c}_{\text{S}}^2/\left(\sqrt{2\pi }{h}_{\text{g}}\right)$ (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 h_d = $h_g\sqrt{\alpha /\left(\alpha +\text{St}\right)}$ (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 = v_α/(St² + l), where smaller particles diffuse more easily than larger ones (Youdin & Lithwick 2007).

To include the effect of dust entrainment in the photoevap-orative wind, we followed the prescription from Sellek et al. (2020), where grains smaller than the entrainment size can be lost with the wind: $$a_{\text{ent}}=\sqrt{\frac{8}{\pi }}\frac{c_{\text{S}}}{G{M}_{*}}\frac{{\dot{M}}_{\text{W}}}{{\Omega }_{*}{\rho }_{\text{s}}},$$(8)

where ${\text{Ω}}_{\ast }=4\pi h_{\text{g}}/\sqrt{h_{\text{g}}^2+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: $${\dot{\text{Σ}}}_{\text{d},\text{w}}(a)=ϵ(a){\dot{\Sigma }}_{\text{g},\text{w}},$$(9)

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).

2.3 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: $$\alpha (r)={\alpha }_0\times \left(1+A_{\text{gap }}\exp \left(-\frac{{\left(r-r_{\text{gap }}\right)}^2}{2w_{\text{gap }}^2}\right)\right),$$(10)

where, α₀ 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) ∝ α₀/α(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 Eq. (7)), though we note that small particles would still be able to pass through the gap by coupling with the gas component.

2.4 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. 2009, 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.

3 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.

3.1 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 Eq. (10)): $${\text{Σ}}_{\text{g}}(r)=\frac{M_{\text{disk }}}{2\pi r_{\text{c}}^2}{\left(\frac{r}{r_{\text{c}}}\right)}^{-1}\exp \left(-r/r_{\text{c}}\right)\frac{{\alpha }_0}{\alpha (r)},$$(11)

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: $$T_{\text{g}}(r)=T_0{\left(\frac{r}{1\text{AU}}\right)}^{-1/2},$$(12)

where T₀ 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 = ϵ₀Σ_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₀ = 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.

3.2 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 pro-toplanetary 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.1 M_* 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² G₀³ (weak), 10³ G₀ (medium), or 10⁴ G₀ (strong; where G₀ 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₀. 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⁷ G₀. Previous work has suggested that the most common experienced UV environment is around 10³ G₀ (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.

3.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_v(a), with v the frequency, we can calculate the optical depth at every radius as: $${\tau }_v={\displaystyle \sum _a{\kappa }_v}(a){\text{Σ}}_{\text{d}}(a),$$(13)

the intensity profile with: $$I_v=B_v(T)\left(1-\exp \left(-{\tau }_v\right)\right),$$(14)

and the total flux as F_v = ∫ I_v dΩ, where dΩ is the differential of the solid angle covered by the disk in the sky, and B_v(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⁻¹ (this is used to calculate the measured disk size r_90%, and the intensity profile I_v).

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: $$T=\sqrt[4]{T_{\text{g}}^4+T_{\text{b}}^4}.$$(15)

4 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.

4.1 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 photoevapora-tion on the gas component is straightforward (Fig. 2, top panel). The gas mass decreases faster for disks undergoing external pho-toevaporation 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 photo-evaporation 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 (Fig. 2, middle panel), and the mass loss in (externally) photoevaporat-ing disks can be described in two stages: a wind dominated loss, and a drift dominated loss (Fig. 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 Myr. 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 photoe-vaporation 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 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 Myr. 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 Myr. 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 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 (Fig. 2, bottom panel) it is interesting to note that this is always below the initial ∈₀ = 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 dust-to-gas ratio than its viscous counterpart, which is due to the reduced 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 Eq. (14)), which are shown in Fig. 5 for the wavelength λ = 1.3 mm. We also indicate the dust masses, the 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 1Myr) we see how the disk subject to external photoevaporation are truncated at r ≈ 60AU–70AU, 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_v (middle panel) follows the trend of the total dust mass. Here, the simulations with photo-evaporative loss show fluxes that range between 0.01–15 mJy after t = 1 Myr, and in particular the simulation with a dust trap manages to display an emission higher than 1 mJy until 4 Myr. 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).

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 rc = 90 AU. The gap substructures are located at rgap = 1/3 rc, and the irradiation field is FUV = 103 G0.

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).

Fig. 4 Evolution of the dust mass loss over time for a photoevaporating disk with FUV = 103 G0. Top: evolution of the dust loss rate 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 Observable quantities derived from the fiducial simulations. Top: intensity profiles 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 Evolution of the disk dust mass for different FUV fields and trap locations (1 /3rc, 2/3rc or none). Top: dust mass remaining in the disk. 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 Evolution of the disk observables, for different FUV fields and trap locations (1 /3rc, 2/3rc or none). Top: evolution disk flux in the 1.3 mm continuum at a distance of 400 pc. Bottom: time evolution of the r90% 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.

4.2 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 (Figs. 6 and 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 was dispersed by the photoevaporative winds within the 5 Myr of simulation time.

For the strongest FUV flux (10⁴ G₀), 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.2Myr until the flux in the millimeter continuum drops.

For the medium FUV flux (10³ G₀), 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 but without external photoevaporation (e.g. Stadler et al. 2022). 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 photoe-vaporative radius reaches the location of the dust trap. This can be seen in Fig. 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² G₀), 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 photoevapo-rative 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 Sect. 4.5.

Fig. 8 Evolution of the disk observables for two different stellar masses, including the gap at rc = 30 AU, with an UV flux of 103 G0. Top: evolution of the flux at the λ = 1.3 mm continuum. Bottom: evolution of the r90% radii (assuming a sensitivity threshold of 0.1 mJy beam-1). The photoevaporative radius is shown for comparison.

4.3 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³ G₀, and fix the dust trap location at 30 AU.

In Fig. 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 occurs 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³G₀. 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 Fig. 9). From this plot, we can expect for disk sizes in pho-toevaporative regions to be determined by their outermost dust trap, independently of the initial extension.

Fig. 9 Same as Fig. 8 for two initial characteristic radius rc. The UV flux is 103 G0.

4.4 Variable UV field

Winter 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.

We conduct a simple simulation with our fiducial parameters (M_* = 1M_*, r_c = 90 AU, r_gap = 30 AU) in which we gradually increase the F_UV irradiation with the following function: $$F_{\text{UV}}=F_{\text{UV},0}+\left(F_{\text{UV},\max }-F_{\text{UV},0}\right){\left(\frac{t}{1\text{Myr}}\right)}^{3/2},$$(16)

where F_UV,0 = 10 G₀ and F_UV,max = 10⁴ G₀ 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 Fig. 10 we show the evolution UV flux and the mass loss rate over time up to 1 Myr.

From the dust distribution shown in Fig. 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-to-centimeter sized pebbles (Fig. 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 (Fig. 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 (Fig. 12), when the flux sharply drops. From our results in Table 3 for the high UV fluxes (10⁴ G₀) 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.

Fig. 10 Time evolution of the UV flux following Eq. (16) (top panel), and evolution of the corresponding dust mass loss rate (bottom panel) by drift (blue) and photoevaporative winds (red) in a disk with variable UV flux. The stellar mass and initial size correspond to the fiducial model (M* = 1 M⊙, rc = 90 AU).

4.5 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: $${\text{Σ}}_{\text{d},\text{ lost }}(r,a)={\displaystyle {\int }_0^{1\text{Myr}}{\dot{\text{Σ}}}_{\text{d},\text{w}}}(r,a)\text{d}t,$$(17)

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² G₀) only grains smaller than 10 µm can be removed from regions beyond 100 AU, while for the strongest irradiation field (10⁴ G₀) 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 Fig. 14 the cumulative size distribution for the parameter space shown in Sect. 4.2 $$M_{\text{dust, lost}}(a)={\displaystyle \sum }^a{\displaystyle {\int }_{r_{\text{in }}}^{r_{\text{out }}}2}\pi r{\text{Σ}}_{\text{d,lost }}(r)\text{d}r,$$(18)

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 (especially in the sub-micron size range). For the weak irradiation environment (10² G₀) 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⁴ G₀) 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), leading to a negative feedback loop where photoevaporation regulates itself. As a first step to study this process we calculate what would be the self-consistent cross section for FUV photons per gas molecule with: $${\sigma }_{\text{FUV}}={\displaystyle \sum _a{ϵ}_w}(a){\kappa }_{\text{UV}}(a)\mu m_{\text{p}},$$(19)

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 Fig. 15 we show the dust-to-gas 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⁻² for all the irradiation regimes and trap configurations. Once the dust has had a chance to grow, the dust-to-gas 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⁻²² cm². This is calculated assuming $${ϵ}_w(a)=\frac{{\dot{\text{Σ}}}_{\text{d}}(a)}{{\dot{\text{Σ}}}_{\text{g}}}.$$(20)

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⁻²² cm²), and find that all our simulations fall below this value, especially 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 Sect. 5.2, we compare the values found for the FUV cross section with previous studies, and discuss if self-shielding by the grains entrained in the wind could be important for the global disk evolution and dispersal process.

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

Fig. 12 Same as Fig. 8, comparing the simulation with constant UV Flux of 103 G0 and the one with increasing UV flux (Eq. (16)).

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 FUV environments, at t = 1 Myr. The disk’s parameters are a star of 1 M⊙, initial size of rc = 90 AU, and a gap substructure at 30 AU. The location of the gap rgap is indicated with a dotted white line.

Fig. 14 Cumulative size distribution of the dust grains lost up to a time of 1 Myr for the simulations shown in Sect. 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 (102 G0) the dust traps have no effect on the lost dust distribution and the three lines overlap.

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

5 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 Sect. 4.1), with depletion timescales that are on the order of 0.1 Myr. 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, especially 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² G₀) to more than 5 Myr, and that in medium radiation environments (10³ G₀) it can increase the disk lifetime to approx. 2 Myr, 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⁴G₀) 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 Sect. 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 HD 147889, and shows signatures of being heated by its bright neighbour.

5.2 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⁻⁴ and a FUV cross section of σ_FUV = 2.7 × 10⁻²³ cm². In Fig. 15 we recalculated the FUV cross section using the dust-to-gas ratio measured directly from the simulation lost material, and found them to be higher (σ_FUV ≈ 10⁻²² cm⁻²) 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, especially 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.

6 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² G₀, the presence of dust traps does prevent drift of dust grains into the star and allow the disk to survive for several Myr. In irradiation environments of 10³ G₀, 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 Sect. 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⁴ G₀ 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, especially 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⁻²² cm² 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.

Acknowledgements

M.G. and P.P. acknowledge funding from the Alexander von Humboldt Foundation in the framework of the Sofja Kovalevskaja Award endowed by the Federal Ministry of Education and Research. P.P. acknowledges funding from UKRI under the UK government’s Horizon Europe funding guarantee from ERC. T.J.H. is funded by a Royal Society Dorothy Hodgkin Fellowship and UKRI guarantee funding for an ERC Consolidator Grant(EP/Y024710/1). S.F. is funded by the European Union under the European Union’s Horizon Europe Research & Innovation Programme 101076613 (UNVEIL). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

References

All Tables

All Figures

Fig. 1 Example mass loss rate grid for a disk orbiting a 1 M⊙ star subject to an irradiation of FUV = 103 G0. The grid determines the gas mass loss ${\dot{M}}_{\text{w}}$ as a function of the gas surface density Σg and radii r, using the FRIED grid from Haworth et al. (2018) and the implementation of Sellek et al. (2020). The solid white line, shows an example gas surface density profile, the dashed vertical line indicates the photoevaporative truncation radius rphot, 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.

In the text

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 rc = 90 AU. The gap substructures are located at rgap = 1/3 rc, and the irradiation field is FUV = 103 G0.

In the text

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).

In the text

Fig. 4 Evolution of the dust mass loss over time for a photoevaporating disk with FUV = 103 G0. Top: evolution of the dust loss rate 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.

In the text

Fig. 5 Observable quantities derived from the fiducial simulations. Top: intensity profiles 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.

In the text

	Fig. 6 Evolution of the disk dust mass for different FUV fields and trap locations (1 /3rc, 2/3rc or none). Top: dust mass remaining in the disk. 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).
In the text

Fig. 7 Evolution of the disk observables, for different FUV fields and trap locations (1 /3rc, 2/3rc or none). Top: evolution disk flux in the 1.3 mm continuum at a distance of 400 pc. Bottom: time evolution of the r90% 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.

In the text

	Fig. 8 Evolution of the disk observables for two different stellar masses, including the gap at rc = 30 AU, with an UV flux of 103 G0. Top: evolution of the flux at the λ = 1.3 mm continuum. Bottom: evolution of the r90% radii (assuming a sensitivity threshold of 0.1 mJy beam-1). The photoevaporative radius is shown for comparison.
In the text

	Fig. 9 Same as Fig. 8 for two initial characteristic radius rc. The UV flux is 103 G0.
In the text

	Fig. 10 Time evolution of the UV flux following Eq. (16) (top panel), and evolution of the corresponding dust mass loss rate (bottom panel) by drift (blue) and photoevaporative winds (red) in a disk with variable UV flux. The stellar mass and initial size correspond to the fiducial model (M* = 1 M⊙, rc = 90 AU).
In the text

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

	Fig. 12 Same as Fig. 8, comparing the simulation with constant UV Flux of 103 G0 and the one with increasing UV flux (Eq. (16)).
In the text

	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 FUV environments, at t = 1 Myr. The disk’s parameters are a star of 1 M⊙, initial size of rc = 90 AU, and a gap substructure at 30 AU. The location of the gap rgap is indicated with a dotted white line.
In the text

Fig. 14 Cumulative size distribution of the dust grains lost up to a time of 1 Myr for the simulations shown in Sect. 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 (102 G0) the dust traps have no effect on the lost dust distribution and the three lines overlap.

In the text

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