© The Authors 2023

1 Introduction

Observations of nearby star-forming regions reveal that the fraction of protoplanetary disks around young stellar objects decreases rapidly with age, indicating that the process of disk dispersal is relatively fast compared to the disk lifetime (Koepferl et al. 2013; Ribas et al. 2015). Theoretical models of disk evolution suggest that photoevaporation could explain this fast dispersal, which occurs when high energy photons in the far ultra-violet (FUV), extreme ultra-violet (EUV), and X-ray wavelength ranges of the spectra hit the gas particles on the disk surface and unbind them from the stellar gravitational potential. When the mass loss rate due to photoevaporation exceeds the local accretion rate, a cavity opens and the disk enters into the photoevaporative dispersal regime, which can clear the remaining material on timescales of 10⁵ yrs from the inside out (e.g., Clarke et al. 2001; Alexander et al. 2006a; Gorti & Hollenbach 2009; Owen et al. 2010; Ercolano & Pascucci 2017).

Photoevaporation has also been proposed as an explanation for transition disks, since the deficit of near-and mid-infrared (NIR and MIR) emission observed in the spectral energy distribution (SED) of these objects is linked to a lack of small grains in the inner regions (Strom et al. 1989; Skrutskie et al. 1990; Espaillat et al. 2014; van der Marel et al. 2016b), which is consistent with a cavity in the inner disk such as those carved by photoevaporative dispersal (Alexander et al. 2006b; Owen et al. 2011; Ercolano et al. 2018), even though not all SED-selected disks might show a cavity in the millimeter continuum, or they might be related to highly inclined disks (van der Marel et al. 2022). Mixed models of photoevaporation and dead zones (Morishima 2012; Gárate et al. 2021) have also succeeded in explaining the high accretion rates found in transition disks (e.g., Cieza et al. 2012; Alcalá et al. 2014; Manara et al. 2014, 2016a,b, 2017), and the presence of a compact inner disk inside their cavities (e.g., Kluska et al. 2018; Pinilla et al. 2019, 2021), which was one of the main limitations of models of photoevaporation acting alone.

However, in addition to the infrared deficits, wide cavities, and high accretion rates, transition disks also seem to be distributed across two populations in terms of their millimeter flux, with a break at approximately F_mm = 30 mJy (d/140 pc)⁻² (Owen et al. 2012; Owen 2016). The fluxes of the population of faint transition disks can be easily reproduced by standard photoevaporation models with dust evolution (Owen & Kollmeier 2019); Gárate et al. 2021), but it is not yet clear if disks undergoing photoevaporation could also reproduce the fluxes of the millimeter bright population.

Observations of rings and gaps in protoplanetary disks (e.g., ALMA Partnership 2015; Andrews et al. 2018; Long et al. 2018; Cieza et al. 2021) indicate that disk density profiles are rich in substructures that can act as dust traps (Whipple 1972; Weidenschilling 1977), and they greatly affect the resulting dust distribution and flux in the millimeter continuum (e.g., Pinilla et al. 2012a,b, 2020). However, models that include consistent dust evolution (i.e., growth, fragmentation, and multiple species, e.g., Birnstiel et al. 2010, 2012; Drążkowska et al. 2019), early substructures¹, along with photoevaporative dispersal have not been widely studied. To our knowledge, only the work of Booth & Owen (2020) has simultaneously considered all three of the mentioned ingredients, in the specific context of the Solar System formation, where the authors show that large amounts of dust can be trapped by Jupiter, therefore decreasing the amount of refractory material delivered to the Sun by the time that photoevaporative dispersal starts to clear out the disk. Thus, it is necessary to further determine the emission in the millimeter continuum of photoevaporating disks where early dust traps are included, and compare them with the fluxes and morphology found in the bright transition disk population.

We note that a common issue of photoevaporative disk models is that these tend to overpredict the fraction of non-accreting transition disks, colloquially dubbed as “relic disks” (Owen et al. 2011), which have not yet been detected by observations. In fact, based on the current observational thresholds, the fraction of relic disks should only be around 3% of the total transition disk population (Hardy et al. 2015), though this fraction should be revised, since several of the observed systems have been identified as non-cluster members in recent years (see Michel et al. 2021). Mechanisms that remove dust grains from the disk, such as radiation pressure (Owen & Kollmeier 2019) or wind entrainment (Franz et al. 2020), could in principle reduce the infrared signal of these relics below detection limits, though it is an open question whether or not these processes would be able to remove enough solid material in a disk where early dust traps were present. Alternatively, it could also be that the fraction of relic disks is simply lower than previously predicted (see Ercolano et al. 2018; Gárate et al. 2021).

In this paper, we studied the evolution of a photoevaporating disk from the point of view of the dust dynamics, using a 1D model. In particular, we focused on how the presence of early substructures (such as the ones caused by planets) affects the resulting dust density and size distribution during photoevapo-rative dispersal, along with the predicted flux in the millimeter continuum (λ = 1.3 mm), and SED in the infrared. In our model we included the growth and fragmentation of multiple dust species (Birnstiel et al. 2010), state-of-the-art models of X-ray and EUV photoevaporation (Picogna et al. 2019), and the loss of dust particles with the photoevaporative winds (Franz et al. 2020).

We further discuss whether the predicted observational signatures of photoevaporating disks can be linked to those observed in transition disks, in terms of millimeter flux and morphology, how the presence or absence of substructures and their properties (location and amplitude) affects the observable features of dispersing disks; and finally, if the dust loss by wind entrainment during the dispersal process can explain why relic disks have not yet been detected.

In Sect. 2, we introduce our disk evolution model and its implementation. In Sect. 3, we present our simulation setup, and the explored parameter space. Our results are shown in Sect. 4, and in Sect. 5 we discuss them in the context of observations. We summarize our results in Sect. 6.

2 Disk model

In this section we present our disk evolution model in 1D, which includes gas and dust advection, dust diffusion, X-ray photoe-vaporation, a prescription for a gap-like substructure, and the evolution of the dust size distribution through coagulation and fragmentation; all assuming that the disk is axisymmetric.

2.1 Gas evolution

The gas evolution is governed by the viscous diffusion and the mass loss due to photoevaporation. Then, the evolution of the gas surface density Σ_g can be described through the following diffusion equation (Lüst 1952; Pringle 1981) $$\frac{\partial }{\partial t}{{\displaystyle \text{Σ}}}_{\text{g}}=\frac{3}{r}\frac{\partial }{\partial r}\left(r^{1/2}\frac{\partial }{\partial r}\left(v{{\displaystyle \text{Σ}}}_{\text{g}}{r}^{1/2}\right)\right)-{{\displaystyle \dot{\text{Σ}}}}_w,$$(1)

where r is the radial distance to the central star and ${\dot{\text{Σ}}}_w$ is the mass loss rate, which depends on the X-ray luminosity L_x from the central star (see Sect. 2.3). The gas viscous evolution is characterized by the kinematic viscosity v, which is defined in Pringle (1981) as $$v=\alpha c_s^2{\text{Ω}}_k^{-1},$$(2)

where α is a dimensionless parameter that represents the magnitude of the gas turbulence (Shakura & Sunyaev 1973), $c_s=\sqrt{k_{B}T/{\text{μm}}_p}$, is the isothermal sound speed, with m_p the proton mass, µ = 2.3 the mean molecular weight, k_B the Boltzmann constant, T the gas temperature. ${\text{Ω}}_k=\sqrt{G{M}_{*}/r^3}$, is the Keplerian orbital speed, with G the gravitational constant, and M_* the mass of the central star.

To induce a gap-like substructure in the gas surface density, such as the one that would be created by a planet, we implemented the following radial turbulence profile that includes a Gaussian bump $$\alpha \left(r\right)={\alpha }_0\times \left(1+A_{\text{gap}}\exp \left(-\frac{{\left(r-r_{\text{gap}}\right)}^2}{2w_{\text{gap}}^2}\right)\right),$$(3)

where α₀ is the base value for the turbulence, r_gap is the location of the gap structure, w_gap is the Gaussian standard deviation that controls the gap width, and A_gap is the amplitude of the bump, which in turn controls the depth of the gap in the gas surface density profile. This Gaussian factor was also used by Pinilla et al. (2020), though in their study it was used to create bumps in the surface density, instead of gaps.

Physically, a local increment in the kinematic viscosity creates a region where the gas diffuses faster, which translates into a gap in the surface density profile. For disks that are in steady state, the gas accretion rate is radially constant, and given by ${\dot{M}}_{\text{g}}=3\pi {\text{Σ}}_{\text{g}}v$ (Pringle 1981). The relation between the gas surface density and viscosity can also be applied to disks that are in quasi-steady state (such as the self-similar solution described by Lynden-Bell & Pringle 1974), where variations in the α turbulence profile yield inversely proportional variations in the gas surface density, approximately following Σ_g ∝ α⁻¹ (see examples in Dullemond et al. 2018; Stammler et al. 2019; Pinilla et al. 2020, among others).

For the gas temperature we assumed that the disk is heated passively by the central star, and therefore its temperature profile is related to the stellar temperature T_* and size R_* through $$T={\theta }_{irr}^{1/4}{\left(\frac{r}{R_{*}}\right)}^{-1/2}{T}_{*},$$(4)

where θ_irr = 0.05 is the disk irradiation angle.

2.2 Dust evolution

For the dust evolution model, we followed the work of Birnstiel et al. (2010), which describes the advection, diffusion, coagulation, and fragmentation of multiple dust species, in response to the interaction with the gas component through the aerodynamic drag force.

The corresponding advection-diffusion equation for the dust surface density Σ_d is $$\frac{\partial }{\partial t}\left(r{\text{Σ}}_{\text{d}}\right)+\frac{\partial }{\partial r}\left(r{\text{Σ}}_{\text{d}}{\upsilon }_{\text{d}}\right)-\frac{\partial }{\partial r}\left(r{D}_{\text{d}}{\text{Σ}}_{\text{g}}\frac{\partial }{\partial r}ϵ\right)=-{{\displaystyle \dot{\text{Σ}}}}_{\text{w},\text{d}},$$(5)

where υ_d corresponds to the dust radial velocity, D_d is the dust diffusivity, ϵ = Σ_d/Σ_g, is the local dust-to-gas ratio, and ${\dot{\text{Σ}}}_{\text{w,d}}$ is the dust loss rate due to wind entrainment (Hutchison et al. 2016; Franz et al. 2020, 2022a,b; Hutchison & Clarke 2021; Booth & Clarke 2021). Here we note that Eq. (5) acts on every individual dust species, and that all dust related quantities are defined as functions of the particle size a.

All components of the dust dynamics for a given particle size are determined by their dimensionless stopping time, the Stokes number $$\text{St}=\frac{\pi }{2}\frac{a{\rho }_s}{{\text{Σ}}_g}\cdot \{\begin{array}{ll}1\hfill & {\lambda }_{\text{mfp}}/a\ge 4/9\hfill \\ \frac{4}{9}\frac{a}{{\lambda }_{\text{mfp}}}\hfill & {\lambda }_{\text{mfp}}/a<4/9.\hfill \end{array}$$(6)

This definition distinguishes between the Epstein drag regime, that dominates when the grain sizes are smaller than the mean free path λ_mfp, and the Stokes regime, that occurs when the particles are large or the gas is very dense (for example in the inner disk).

The radial velocity of a dust particle is then $${\upsilon }_{\text{d}}=\frac{1}{1+{\text{St}}^2}{\upsilon }_v-\frac{2\text{St}}{1+{\text{St}}^2}\eta {\upsilon }_k,$$(7)

following Nakagawa et al. (1986) and Takeuchi & Lin (2002), where υ_v is the gas advection velocity due to viscous diffusion and η = − (1/2) (h_g/r)² dlnP/dln r, is the relative difference between the Keplerian velocity υ_k and the gas orbital velocity, due to its own pressure support. The isothermal pressure is defined as $P={\rho }_{\text{g}}{}_{,0}{c}_s^2$, with ρ_g,0 the gas volume density at the midplane, and $h_{\text{g}}=c_s{\text{Ω}}_k^{-1}$, is the gas scale height.

Finally, the dust diffuses with a diffusivity D_d = ν/(1 + St²) (Youdin & Lithwick 2007), where we note that D_d ≈ v for particles with St ≪ 1.

In addition to advection and diffusion, dust species also grow and/or fragment depending on their relative velocities and their collision rate, where the evolution of the grain size distribution is governed by the Smoluchowski equation (Birnstiel et al. 2010). In a typical protoplanetary disk there are two regimes of dust growth: fragmentation limited, or drift limited (Birnstiel et al. 2010, 2012).

The fragmentation limit occurs when the collision velocities of larger dust grains surpass the fragmentation velocity threshold υ_frag (which depends on the material properties), resulting in destructive collisions, and replenishing the population of small grains (e.g., Ormel & Cuzzi 2007; Brauer et al. 2008; Birnstiel et al. 2009). This regime can be typically found in the inner regions of protoplanetary disks, regions with high turbulence, and in pressure maxima, where the drift velocity is zero, and is given by $${\text{St}}_{\text{frag}}=\frac{1}{3}\frac{{\upsilon }_{\text{frag}}^2}{\alpha c_s^2}.$$(8)

The drift limit, on the other hand, occurs when the dust grains grow to the point where the drift timescale is shorter than the growth timescale (Birnstiel et al. 2010). This regime appears in regions with steep pressure gradients, such as the outer disk and regions with low dust-to-gas ratios, where the maximum size that a grain can grow to is approximately $${\text{St}}_{\text{drift}}={\left|\frac{\text{dln\hspace{0.17em}}P}{\text{dln\hspace{0.17em}}r}\right|}^{-1}\frac{{\upsilon }_k^2}{c_s^2}{ϵ}_{\text{tot}},$$(9)

where ϵ_tot refers to the local dust-to-gas ratio of all dust species combined.

2.3 Photoevaporation model

When high energy radiation from the central star hits the disk surface layer, the material is heated and unbound from the stellar gravitational potential in a process called photoevaporation, which ultimately leads to disk dispersal (Clarke et al. 2001; Alexander et al. 2006a,b; Alexander & Armitage 2007).

In our model, we implemented the mass loss rate profile from Picogna et al. (2019, see their Eqs. (2)–(5)) into the sink term ${\dot{\text{Σ}}}_w$ of the gas diffusion equation (Eq. (1)) to simulate the effect of X-ray and EUV photoevaporation, where the total mass loss rate increases with the stellar X-ray luminosity L_x.

Figure 1 shows an example of the ${\dot{\text{Σ}}}_w\left(r\right)$ profile, that distinguishes between the case when the disk is still young and without a cavity, and the case after the photoevaporative cavity opens, where the cavity edge is directly irradiated by the central star. This model is valid for a 0.7 M_⊙ star irradiating a disk with the X-ray spectrum as given by Ercolano et al. (2009). Ercolano et al. (2021) and Picogna et al. (2021) expand this model to a range of stellar masses and apply observationally determined X-ray spectra. Given that the new models are qualitatively similar to those of Picogna et al. (2019) used here, we do not expect that their implementation in our work would lead to significant changes in our conclusions.

From Eq. (1), we see that the gas evolution can be dominated by either viscous diffusion or by photoevaporation. In the early stages, when the disks are more massive, the viscous accretion is generally thought to be the dominant evolution mechanism, while photoevaporation will dominate the disk dispersal in later stages, removing the remaining material by opening a cavity from the inside out (Clarke et al. 2001).

Along with the gas removal through photoevaporative winds, we can also expect for a fraction of the dust grains to be entrained in the photoevaporative flow (Hutchison et al. 2016; Franz et al. 2020). To quantify the dust loss rate we define a sink term in Eq. (5) $${{\displaystyle \dot{\text{Σ}}}}_{\text{w},d}={ϵ}_{\text{w}}{{\displaystyle \dot{\text{Σ}}}}_{\text{w}},$$(10)

where ϵ_w represents the dust-to-gas loss ratio, and is defined in our model as the mass fraction of particles that are small enough to couple to the gas motion with a ≤ a_w, and that lie above the wind launching surface with z ≥ h_w, where h_w and a_w are free parameters in our model.

The final ingredient to find ϵ_w is to define a vertical structure for the gas and dust volume densities ρ_g and ρ_d. Following Fromang & Nelson (2009), and assuming that the gas is in vertical hydrostatic equilibrium (with constant temperature in the vertical direction), we model the vertical density distribution as $${\rho }_{\text{g},\text{d}}\left(z\right)=\frac{{\text{Σ}}_{\text{g},\text{d}}}{\sqrt{2\pi }{h}_{\text{g,d}}}\exp \left(-\frac{z^2}{2h_{\text{g},\text{d}}^2}\right),$$(11)

where the dust scale height h_d is $$h_{\text{d}}=h_{\text{g}}\cdot \min \left(1,\sqrt{\frac{\alpha }{\min \left(\text{St,1/2}\right)\left(1+{\text{St}}^2\right)}}\right),$$(12)

following Youdin & Lithwick (2007) and Birnstiel et al. (2010), where we notice that the dust scale height is smaller than the gas one, since large grains with (St ≳ α) tend to settle toward the midplane.

The formal expression for the dust-to-gas loss ratio can then be written as $${ϵ}_{\text{w}}\left(a\right)=\frac{{\int }_{h_{\text{w}}}^{\infty }{\rho }_{\text{d}}\left(z,a\right)\text{d}z}{{\int }_{h_{\text{w}}}^{\infty }{\rho }_{\text{g}}\left(z\right)\text{d}z},$$(13)

which due to settling and growth is always smaller than the local dust-to-gas mass ratio ((ϵ_w < ϵ).

Our photoevaporation model assumes that all the small grains above the wind launching surface are fully entrained and that none of the removed material (gas or dust) falls back onto the disk (Clarke & Alexander 2016); Picogna et al. 2019; Franz et al. 2020; Sellek et al. 2021).

Fig. 1 Example of the gas loss rate profile, following the photoevaporation model from Picogna et al. (2019, their Eqs. (2) and (4)), for Lx = 1030 erg s−1. The figure shows the mass loss rate before the photoevap-orative cavity opens (solid line) and after the photoevaporative cavity opens (dashed line), with the cavity edge located at rcavity ≈ 30 AU after 1.9 Myr of evolution.

3 Simulation setup

We performed our simulations using the code DustPy²,³ (Stammler & Birnstiel 2022a), that solves the gas and dust surface density evolution (Eq. (1) and Eq. (5)) following the model from Birnstiel et al. (2010).

We implemented the photoevaporation model described in Sect. 2.3 and induced gap-like substructures using the α turbulence profile from Eq. (3). For our study, we focused on the impact that the X-ray luminosity L_x, the gap amplitude A_gap, and the gap location r_gap have on the disk evolution.

For each simulation we tracked the evolution of the disk mass M_g,d (in gas and dust), photoevaporative cavity size r_cavity, flux in the millimeter continuum F_mm at 1.3 mm, and SED. Our goal is then to compare these properties to the values observed in proto-planetary disk populations. In particular, we want to determine if these photoevaporating disks could be bright enough to explain the transition disk millimeter fluxes, or if they are too faint to be detected at all, and therefore related to the relic disk problem (Owen et al. 2011, 2012; Owen & Kollmeier 2019).

In this section we describe the initial conditions and numerical grid setup, the radiative transfer model that we used to post-process the simulations and obtain the millimeter fluxes and SEDS, and finally our parameter space exploration, which we used to study the impact of the X-ray luminosity and gap properties.

3.1 Initial conditions and numerical grid

In our simulations, the central star has a mass of M_* = 0.7 M_⊙, a radius of 1.7 R_⊙, and a temperature of 4500 K, selected to match the stellar parameters from the photoevaporative models of Picogna et al. (2019) and Owen & Kollmeier (2019). For these stellar properties, the disk has a temperature of ≈190 K at 1 AU (see Eq. (4)).

For the initial gas surface density profile we used a modified version of the Lynden-Bell & Pringle (1974) self-similar solution $${\text{Σ}}_{\text{g}}\left(r\right)=\frac{M_{\text{disk}}}{2\pi r_c^2}{\left(\frac{r}{r_c}\right)}^{-1}\exp \left(-r/r_c\right)\frac{{\alpha }_0}{{}^{a\left(r\right)}},$$(14)

where M_disk = 0.05 M_⊙ is the total disk mass, and r_c = 60 AU is the disk characteristic radius.

To introduce a gap in the surface density profile from the beginning of the simulation, we added the factor α0/α(r) to the self-similar solution. Then, the resulting gap structure is consistent and sustained by the turbulence profile defined in Eq. (3) (Stammler et al. 2019; Pinilla et al. 2020; Stadler et al. 2022), where the gap center is located at r_gap and the amplitude (i.e., the depth of the perturbation in the surface density profile) is determined by A_gap. In the case where the gap amplitude is A_gap = 0 (i.e., no gap), we recover the traditional self-similar solution. We used a value of α₀ = 10⁻³ for the disk base turbulence.

The initial dust-to-gas ratio is ϵ₀ = 1.5 × 10⁻², and the initial dust size distribution follows the MRN distribution (Mathis et al. 1977), with an initial maximum grain size of a₀ = 1 µm. Our initial dust-to-gas ratio is higher than the canonical 1%, motivated by a recent study by Lebreuilly et al. (2020), which indicates that protoplanetary disks may heritage higher dust-to-gas ratio than the ISM from the protostellar collapse.

For the dust grains we assumed that these are compact and covered by ice, with a material density of ρ_s = 1.6 g cm⁻³ and a fragmentation velocity of υ_frag = 10 m s⁻¹ (Wada et al. 2011; Gundlach et al. 2011; Gundlach & Blum 2015), though notice that recent results suggest that the fragmentation velocity of icy grains could be lower than previously thought (Gundlach et al. 2018; Musiolik & Wurm 2019; Steinpilz et al. 2019).

We used a logarithmically spaced radial grid going from 4 AU to 300 AU with n_r = 200 radial cells, and a logarithmically spaced mass grid from 10⁻¹² g to 10⁵ g (approx. 0.5 µm to 20 cm in grain sizes) with n_m = 120 cells. Finally, in order to determine the fraction of dust that is entrained in the photoevaporative wind ϵ_w, we employed a 1+1D approach in which we constructed a vertical grid locally at every radial grid cell to solve the integrals in Eq. (13). This grid is defined in function of the gas scale height, going from the midplane to 10 h_g with n_z = 100 cells.

We saved the simulation outputs every 0.1 Myr, and terminated the simulations when the photoevaporative cavity exceeded 120 AU in size, since other photoevaporation regimes are more likely to become dominant over the X-ray driven dispersal for larger cavity sizes (e.g. FUV; Gorti & Hollenbach 2009).

3.2 Radiative transfer and optically thin approximation

To obtain the millimeter fluxes F_mm at λ = 1.3 mm we can take two approaches: the vertical slab approximation, or the complete radiative transfer calculation. For the vertical slab approximation we used the vertically integrated surface density to calculate the optical depth τ_v = Σ_a κ_v(α)Σ_d(a), where κ_v(a) is the absorption opacity and v is the frequency, and obtain the total flux at λ = 1.3 mm (v = 230 GHz) with $$F_{\text{mm}}={\displaystyle \int B_v\left(T\right)\left(1-\exp \left(-{\tau }_v\right)\right)\text{dΩ,}}$$(15)

where B_v is the Planck function, dΩ is the solid angle differential, and T is the vertically isothermal dust temperature (Eq. (4)). This approach is ideal to quickly compute the fluxes for all snapshots directly from the dust distribution, but has the drawback that it is only reliable for low optical depths, it neglects the effect of self-scattering, and that the temperature profile may not be consistent with that of a irradiated disk, especially at the edge of the photoevaporative cavity.

To obtain more accurate fluxes for key snapshots, we used the radiative transfer code RADMC-3D⁴ (Dullemond et al. 2012), to recalculate the dust temperature, the millimeter fluxes at 1.3 mm, and the SED between 0.1 µm and 1 cm. For our calculations we considered the complete treatment of scattering, that includes polarization and anisotropy (Kataoka et al. 2015). The radiative transfer is performed on a azimuthally symmetrical spherical grid, where the radial coordinate matches the logarithmically spaced grid from Dustpy (see Sect. 3.1), and the colatitude coordinate covers the entire domain with n_θ = 180 grid cells. We used 10⁷ photon packages to calculate the thermal structure, 2.5 × 10⁵ photon packages to calculate the emission in the millimeter continuum, and 10⁴ photon packages to calculate the SEDs. To account for the full scattering treatment we subdivided the azimuthal coordinate in n_ϕ = 64 grid cells. For this work we also assumed that our disks are face-on (the inclination is i = 0), and that they are located at a distance of d = 140 pc, which is the typical distance of the nearby star-forming regions (Dzib et al. 2018; Roccatagliata et al. 2020).

For both the optically thin and the radiative transfer setup we used the opacity model from the DSHARP survey (Birnstiel et al. 2018), which assumes compact grains composed of water ice, troilite, refractory organics, and astronomical silicates (Henning & Stognienko 1996); Draine 2003; Warren & Brandt 2008). Then we used the code OpTool⁵ (Dominik et al. 2021) to obtain the opacities, following the Mie theory for compact grains, for all 120 grain sizes tracked by the Dustpy simulations.

While the DSHARP opacities provide a convenient framework that is common to several recent studies, it is not clear whether these represent the true absorption of dust grains accurately. For example, the model from Ricci et al. (2010), based on the optical constants of Zubko et al. (1996), Draine (2003), and Warren & Brandt (2008), leads to absorption opacities that are approximately one order of magnitude higher than the DSHARP opacities in the millimeter continuum, which leads to higher optical depths and fluxes (Zormpas et al. 2022; Stadler et al. 2022). Since the fluxes obtained from radiative transfer calculations are dependent on the selected opacity model, we included in our results a comparison between the fluxes obtained with the Birnstiel et al. (2018) and the Ricci et al. (2010) opacity models.

3.3 Parameter space: X-ray luminosity and gap properties

In this study we want to understand what effect the presence or absence of substructure has on the disk observable quantities during photoevaporative dispersal. To explore the parameter space we selected two fiducial simulations: one without a gap (A_gap = 0), and one with a gap (A_gap = 4, i.e., a decrease by a factor of 0.2 in the local Σ_g) located at r_gap = 40 AU, where both simulations have the X-ray luminosity L_x = 10³⁰ erg s⁻¹. For reference, a gap located at 40 AU and with an amplitude of 4, is what we would expect from a planet of 225 M_⊕ (approx. twice the mass of Saturn), or planet-to-star mass ratio of q ≈ 9.5 × 10⁻⁴, following the Kanagawa et al. (2017) gap model.

Afterward, we repeated our study for different X-ray luminosities, while keeping the fiducial gap properties. Finally we studied the effect of the different gap locations and amplitudes, this time keeping the fiducial X-ray luminosity. Table 1 shows the X-ray luminosity and gap properties of our parameter space, with the fiducial values in boldface.

The X-ray luminosities were selected from within the range of the Taurus luminosity distribution (Preibisch et al. 2005). Each value of L_x can also be understood in terms of the resulting total mass loss rates ${\dot{M}}_w$, which are respectively 4.6 × 10⁻⁹, 1.6 × 10⁻⁸ and 3.2 × 10⁻⁸ M_⊙ yr⁻¹ for the values listed in Table 1 (see Eq. (5) from Picogna et al. 2019).

The gap locations were selected to be within (or at) the disk characteristic size r_c = 60 AU, and the maximum gap amplitude was selected to ensure that dust trapping is effective at the different gap locations. For the gap widths, we chose to apply a simple prescription of $$w_{\text{gap}}=5\left(\frac{r_{\text{gap}}}{40\text{AU}}\right)\text{AU,}$$(16)

which roughly matches the widths of the dust traps from Pinilla et al. (2020), and is always larger than the local scale height for our parameter space.

Finally, the amount of dust removed by photoevaporation in our model depends both on the maximum entrainment size and on the scale height of the wind launching region. As fiducial values for our simulations, we assumed that only particles smaller than a ≤ a_w = 10 µm (Hutchison et al. 2016; Franz et al. 2020; Hutchison & Clarke 2021; Booth & Clarke 2021) can be carried by the photoevaporative winds, and that the particles must be at least above z ≥ h_w = 3h_g, though the photoevaporative surface can be located at higher altitudes.

In Appendix A we further explored the parameter space for aw and h_w, though we do not expect for the resulting dust distribution to be greatly affected by the exact parameter values, since due to grain growth and settling dust loss should be small in comparison to the gas loss (i.e., ϵ_w < ϵ; Franz et al. 2022a). Motivated by the possibility of more efficient dust removal due to additional mechanisms such as radiation pressure (Owen & Kollmeier 2019), an FUV component in the wind (e.g., Gorti & Hollenbach 2009), or grains lifted by magneto-hydrodyamical (MHD) winds (Miyake et al. 2016), we also included a model in which all the dust is fully entrained in the wind (ϵ_w = ϵ). We also included a comparison of the dust mass evolution of our fiducial model against the entrainment prescription of Booth & Clarke (2021).

4 Results

4.1 Fiducial models

In this section we present our results for the evolution of photo-evaporating disks, focusing on the effect that early substructures (represented through a primordial gap in the gas component) have on the dust component, in terms of the distribution of solids and the corresponding observable quantities (millimeter fluxes and SEDs).

We often refer to the disk without the primordial gap as a “smooth disk”, and the disk with the primordial gap as a “structured disk”. We also distinguish between the “gap” structure that is created through the variation in the α viscosity profile (Eq. (3)), and the “cavity” that is carved by photoevaporative dispersal.

4.1.1 Evolution of the dust distribution

A disk with a gap-like substructure can efficiently trap dust grains at the local pressure maximum, so long as the substructure forms before any significant radial drift occurs (Stadler et al. 2022). In contrast, in a disk without substructure the dust drifts very efficiently toward the star (Birnstiel et al. 2010; Pinilla et al. 2012b, 2020). Our fiducial models show that, by the time that photoevaporation starts clearing the gas component from the inside out, after ~ 1 Myr of disk evolution (for L_x = 10³⁰ erg s⁻¹), the structured disk retains a higher mass of solids than the smooth disk, 77 M_⊕ and 4 M_⊕ respectively, out of the initial ~ 200 M_⊕ (see Fig. 2, top panel).

Once photoevaporation opens up a cavity in the inner regions, further dust drift toward the star is completely halted, since the edge of the cavity is also a pressure maximum that can trap solids (see the evolution of the surface densities, Fig. 3). From this point onward, the dust loss is driven exclusively by the entrainment in the photoevaporative winds, though additional loss terms such as planetesimal formation and removal by radiation pressure are neglected in our model (Stammler et al. 2019; Owen & Kollmeier 2019, see discussion in Sect. 5.3).

Figure 2 (bottom panel) shows that the dust loss rate by pho-toevaporative entrainment increases between one and two orders of magnitude after the photoevaporative cavity opens. One reason is that the gas loss rate locally increases at the edge of the photoevaporative cavity, where the material is directly irradiated by the central star. In Fig. 1 we see how the gas loss profile changes in this “open cavity” scenario, with a sharp spike at the location of the cavity edge (see also Picogna et al. 2019, their Eq. (4)). The other reason, and perhaps more importantly, is that dust growth becomes limited by fragmentation around the pressure maxima which replenishes the population of small grains that are more easily entrained with the wind, and the gas surface density is also reduced at the photoevaporative cavity edge, sharply reducing the maximum grain size. In contrast, the growth in the gas rich regions with steeper pressure gradients is limited only by drift, which results in a dust distribution dominated by large particles (this can bee seen from the grain size distributions in Figs. 4 and 5), that are not easily entrained with the photoevaporative wind.

After the photoevaporative cavity opens, the remaining dust mass decreases from 77 M_⊕ to 55 M_⊕ for the structured disk model, and from 4 M_⊕ to 3 M_⊕ for the smooth disk model. These values imply that the total dust loss across the disk lifetime (or at least until the cavity size reaches r_cavity = 120 AU in our model) is mostly dominated by drift during the early stages of disk evolution, rather than entrainment in the photoevaporative winds (see also Ercolano et al. 2017).

From the surface density profiles (Fig. 3) and the grain size distributions (Figs. 4 and 5), we find that once a photoevaporative cavity opens, the remaining solid material is dragged along with the cavity outer edge, following the moving pressure maximum. For the smooth disk, this leads to the formation of a single dust trap at 1.4 Myr, that moves outward as time passes. For the structured disk, on the other hand, we find that between 1.2 and 1.4 Myr there are two traps present, one that follows the photoevaporative cavity, and the other at the outer edge of the primordial gap, which should lead to a distinct disk morphology featuring two rings. Eventually, both dust traps merge into one when the cavity catches up with the gap location, which then continues to move outward. We infer that the two ring morphology is more likely to be observed if the primordial dust trap is located at larger radii than the photoevaporative cavity opening radii. The latter would delay the merging of the two rings and increase the window of observation, though to get an accurate estimate of the likelihood of observing this evolutionary stage, a population synthesis model would be required. We also note that if a dead zone is present in the inner disk, the photoevaporative cavity opening radius can be located beyond 10 to 20 AU (Gárate et al. 2021), meaning that disks with primordial dust traps located inside the dead zone radius (such as Jupiter’s current orbit) would not lead to the described two ring morphology.

This particular behavior in the evolution of structured disks leads to a degeneracy between the properties of an observed dust ring and its potential origin, which is of particular interest for the study of transition disks, and we discuss more about it in Sect. 5.2.

Fig. 2 Top: Mass evolution of the disk mass in gas (red lines) and dust (blue lines). Bottom: Evolution of the mass loss rate of gas and dust by photoevaporative winds. The markers indicate the moment when photo-evaporation opens a cavity in the inner disk (“+” for the initially smooth disk, “x” for the initially structured disk).

Fig. 3 Evolution of the surface density profiles disks from t = 1 Myr, and plotted every 0.1 Myr (solid lines, with line opacity increasing with time). The dust surface density accounts for all the grain sizes. The initial condition is shown with dashed lines.

Fig. 4 Dust size distribution for a smooth disk at 1.2, 1.6, and 2.1 Myr. The drift (cyan) and fragmentation (pink) growth limits are also indicated.

Fig. 5 Dust size distribution for a structured disk at 1.1, 1.3, and 1.6 Myr. The drift (cyan) and fragmentation (pink) growth limits are also indicated. The last two snapshots were selected to match the ones of Fig. 4 in terms of the photoevaporative cavity size.

4.1.2 Millimeter emission and SED

Because the dust masses between the structured and smooth disks differ by over an order of magnitude during the photoevaporative dispersal, this also results in a similar difference between their corresponding luminosities in the millimeter continuum (Fig. 6). In our models, the flux of the structured disk is F_mm ≈ 65 mJy (obtained from Eq. (15)) by the time the cavity opens, and remains approximately constant until the cavity reaches the location of the primordial dust trap at r ≈ 50 AU, afterward the flux continues to decrease, reaching 44 mJy by the time the cavity has grown to 100 AU. The smooth disk flux simply decreases approximately from 4 mJy to 2 mJy.

We also note that the radiative transfer calculations with RADMC-3D differ only by a small factor from the values obtained with the vertical slab approximation when using the DSHARP opacity model. The difference in both fluxes is likely due the direct heating of the photoevaporative cavity edge by the stellar irradiation, and the proper treatment of the scattering and optical depth.

Finally, the SEDs (Fig. 7) show a similar behavior for both the smooth and structured disk, where the deficit in the NIR to MIR wavelengths becomes more prominent as the cavity size grows, and small grains are removed from the inner regions, meaning that our models would be classified as transition disks by their SED (Espaillat et al. 2014). Additionally, the disks in our model display a high emission in the far infrared (FIR, around 100 µm) which are comparable to those of the transition disk Sz-Cha (van der Marel et al. 2016b; Gaia Collaboration 2020, with the distance rescaled to 140 pc), though we remark that this is not intended to be a representative comparison with the transition disk population. In Sect. 5.3 we discuss the implications of the FIR excess in the context of the relic disk problem.

Fig. 6 Millimeter fluxes Fmn, as a function of the size of the photoevaporative cavity for the smooth (dotted) and structured (solid) disk models. The solid and dotted lines represent flux from the optically thin approximation (Eq. (15)). The markers are the fluxes obtained with RADMC-3D (“+” for the smooth disk, “x” for the structured disk). The disks are assumed to be at 140 pc, and no inclination. The cavity size measurement is based on the dust distribution of millimeter sized grains.

Fig. 7 SEDs for the smooth (dotted) and structured (solid) disk models, when the photoevaporative cavity size is rcavity ≈ 15 AU (gray), and rcavity ≈ 100 AU (black), assuming a distance of 140 pc and face-on. The data points show the SED of SzCha from van der Marel et al. (2016b) survey, re-scaled to a distance of 140 pc.

Fig. 8 Evolution of the dust mass (top) and disk flux at λ = 1.3 mm (bottom, assuming a distance of 140 pc), for different X-ray luminosities Lx, with the black line corresponding to the fiducial value. The markers indicate the moment when photoevaporation opens a cavity in the inner disk (“+” for the smooth disk, “x” for the structured disk).

4.2 Effect of the X-ray luminosity

In this section we test the impact that different X-ray luminosities have on the evolution of the disk mass in the dust component, and on the corresponding flux in the millimeter continuum (see Fig. 8). For the structured disk, we use again the fiducial gap amplitude of A_gap = 4 located at r_gap = 40 AU.

We notice that for higher X-ray luminosities the photoevaporative cavity opens earlier (due to the higher mass loss rates), and that both the dust mass and millimeter flux are also higher when the photoevaporative cavity opens. This occurs because the dust drift, dust diffusion, and wind entrainment processes had less time to remove solid material from the disk.

The difference can be clearly seen in the mass and flux evolution of the smooth disk, where dust drift can only be stopped by the pressure maximum corresponding to the photoevaporative cavity. For the highest X-ray luminosity (L_x = 3 × 10³⁰ erg s⁻¹) the millimeter flux of the smooth disk is F_mm = l2 mJy (at 0.8 Myr, when the cavity opens), while for the lowest X-ray luminosity (L_x = 3 × 10²⁹ erg s⁻¹) the flux is only F_mm = 0.4 mJy (at 4.3 Myr).

On the other hand, while the structured disk can trap dust particles at the local pressure maximum, some of them will still diffuse through the gap and be lost to the star. We find that dust entrainment only accounts for a minor fraction of the dust removal in structured disks before the opening of the photoevaporative cavity, with rates between 10⁻⁶ M_⊕ yr⁻¹ and 10⁻⁵ M_⊕ yr⁻¹ (for L_x = 3 × 10²⁹ erg s⁻¹ and 3 × 10³⁰ erg s⁻¹, respectively). The millimeter flux of the structured disk is F_mm = 70 mJy when the cavity opens (at 0.7 Myr) for the highest X-ray luminosity, and F_mm = 49 mJy for the lowest X-ray luminosity (at 3.6 Myr). We expect that disks with multiple dust traps would be able to retain more material and prevent further dust loss due to diffusion (Pinilla et al. 2012b, 2020).

Another feature that we observe for each pair of simulations with the same X-ray luminosity, is that the inner cavity opens earlier in the structured disk than in the smooth one by 0.2 Myr, a behavior that is also seen in the simulations of Rosotti et al. (2013) when the effects of photoevaporation and planet-disk interactions are considered. This occurs because the presence of the gap structure seems to speed up the viscous evolution of the disk by a small factor, reducing the gas accretion rate in the inner regions faster, and allowing for photoevaporative dispersal to start earlier. This is also the reason why the gas mass decreases slightly faster in the initially structured disk than in the initially smooth disk shown in Fig. 2 (top panel).

Fig. 9 Same as Fig. 8, but for structured disks with different gap amplitudes, with a fixed location at rgap = 40 AU and X-ray luminosity of Lx = 1030 erg s−1. Notice that the axis scales are different from Fig. 8.

4.3 Effect of the trap location and amplitude

In this section we test the impact that the gap amplitude and location have on the evolution of the mass in the dust component and the respective flux in the millimeter continuum (Figs. 9 and 10). For the X-ray luminosity we use the fiducial value of L_x = 10³⁰erg s⁻¹.

We observe that a minimum amplitude of A_gap = 2, which roughly corresponding to a Saturn mass in our models, seems to be required for a gap at r_gap = 40 AU to effectively trap the dust particles that drift from the outer region (Fig. 9, top panel). The dust mass for this disk, when the inner cavity opens, is about 60 M_⊕, in comparison with the 77 M_⊕ measured for the disk with the deeper gap (A_gap = 4). In contrast, the disk with A_gap = 1 is unable to stop the dust drift and is almost indistinguishable from the completely smooth disk, retaining a dust mass of only 7 M_⊕ at the cavity opening.

The location of the gap (relative to the disk characteristic radius r_c) on the other hand determines the size of the dust reservoir that can be retained at a given dust trap. Because dust grains tend to drift inward, all the dust that is initially inside the gap location is rapidly lost into the star, as can be seen during the first 0.1 Myr of disk evolution (Fig. 10, top panel). Only the dust grains that are further away than the local pressure maximum can be potentially trapped in it (if we ignore diffusion and wind entrainment), and therefore gaps that are closer to the star result in a higher dust content during the disk evolution and dispersal. For comparison, the dust mass of the disk with the innermost gap at 20 AU, is approximately double the mass of the disk with the outermost gap at 60 AU, at the moment of the cavity opening.

In terms of the millimeter flux we notice that the disk with r_gap = 20 AU is slightly fainter than the fiducial model with r_gap = 40 AU upon the opening of the inner cavity. This might sound counter-intuitive, since the simulation with r_gap = 20 AU has a larger dust mass (Fig. 10) at all times. However, we also notice that a dust trap located further inside into the disk will have both a higher optical depth (τ_v > 1) and a smaller surface area, where both of these effects would then contribute to reduce the total flux observed (see Eq. (15)). In other words, at the moment of the cavity opening the disk with the innermost dust trap has more solid material, but this material is “hidden” from the observer. As the cavity starts growing and the dust trap moves out, the disk with the initially innermost gap becomes the brightest, reaching F_mm = 78 mJy at its peak (after the inner cavity opens).

We note that a disk with multiple traps located in the inner and outer regions should be able to retain a higher fraction of the initial dust mass, and display higher fluxes in the millimeter continuum. We also want to point out that the ratio between the gap location and the disk initial size (r_gap/r_c) should be more important to predict the amount of dust trapped than their absolute value, since that the same values of r_gap presented in this section may trap more (or less) solid material for disks that are initially more extended (or compact).

Fig. 10 Same as Fig. 8, but for structured disks with different gap locations, with a fixed amplitude of Agap = 4 and X-ray luminosity of Lx = 1030 erg s−1. Notice that the y-axis scale is different from Fig. 9.

4.4 Comparison between opacities

Due to the uncertainty in the opacities of dust grains in pro-toplanetary disks, particularly in the amount of carbonaceous material, we chose to perform an additional set of radiative transfer calculations, this time with the opacity model of Ricci et al. (2010), and compare the corresponding millimeter fluxes against those obtained with the Birnstiel et al. (2010; DSHARP) model in Sect. 4.1.2.

From Table 2 we observe that the flux in the millimeter continuum is always higher when the Ricci et al. (2010) opacity model is considered, which is to be expected since the absorption opacities are also higher. The trend between smooth and structured disks is maintained no matter the opacity model used, with the smooth disks being fainter (F_mm < 20 mJy) and the structured disks being brighter (F_mm > 110 mJy).

Finally, we notice that the relative increase in flux when using the Ricci et al. (2010) opacities is higher in the smooth disks (a factor between 4 to 5) than in the structured disks (a factor between 2 to 4). This occurs because the structured disks have regions that are optically thick (at the dust traps, for example), while on the other hand most of the smooth disks are optically thin, and therefore more sensitive to changes in the opacity model.

5 Discussion

5.1 Explaining the observed transition disks with photoevaporation models

Since first discovered, transition disk properties have remained a challenge to theoretical models. These objects display deep cavities in the dust component, as probed by the deficit in NIR and MIR emission (Strom et al. 1989; Skrutskie et al. 1990; van der Marel et al. 2016b) and resolved continuum observations (Andrews et al. 2011; Pinilla et al. 2018; Francis & van der Marel 2020), while displaying relatively high accretion rates (${\dot{M}}_{\text{acc}}\sim {10}^{-10}$ to 10⁻⁸ M_⊙ yr⁻¹) that indicate the presence of long lived inner disks (Cieza et al. 2012; Alcalá et al. 2014; Manara et al. 2017). Additionally, transition disks seem to be distributed across a mm-faint (F_mm < 30 mJy) and a mm-bright (F_mm > 30 mJy) population, according to their fluxes in the millimeter continuum (see the review by Owen 2016). However, it should be noted that the faint transition disks, which are identified based on their SEDs, could be instead highly inclined disks rather than actual transition disks with cavities (van der Marel et al. 2022).

Giant planets have often been proposed as an explanation for transition disk properties, since these can carve deep gaps in the millimeter continuum by trapping millimeter grains, while still allowing for the flow of gas toward the inner disk (Dong et al. 2012; Pinilla et al. 2012a; de Juan Ovelar et al. 2013; Owen 2014). On that line, a recent study of van der Marel & Mulders (2021) proposes that the population of observed transition disks can be linked to the population of detected exoplanets, under the assumption that planets more massive than Jupiter are responsible for the transition disk substructures, and that these planets then migrate inward.

Yet, despite the high masses predicted for the planetary companions (Jupiter-sized or larger, Muley et al. 2019; van der Marel et al. 2021), to this day the only disk with a confirmed planet detection is PDS-70 (Keppler et al. 2018). While some of these hypothetical planets would still be beyond the current detection limits (Asensio-Torres et al. 2021), we cannot discard that other mechanisms may be actually responsible for transition disks.

The simulations presented in this work (and in combination with Gárate et al. 2021), offer an additional pathway to explain the properties of transition disks without the presence of very massive planets, which would also explain why we have not detected these giants in the first place. In Sect. 4 we showed that when substructures are present, the dust grains are retained in the local pressure maxima through the lifetime of the protoplanetary disk. Then, once the inside-out dispersal due to photoevaporation begins, all the trapped material is dragged along with the edge of the photoevaporative cavity, which results in a bright disk with a wide expanding cavity (Fig. 6), consistent with the mm-bright population of transition disks.

Now, notice that the only requirements imposed on these substructures to explain transition disk properties of the mm-bright population would be a minimum amplitude (i.e., gap depth) to trap solid material, and dust traps located closer to the star favor would favor brighter disks (see Figs. 9 and 10). This means that even planets with masses between those of Saturn and Jupiter, without imposing strong constraints on their location, could create the necessary substructures to reproduce a population of millimeter bright transition disk, under the assumption that internal photoevaporation is the primary dispersal process responsible for opening the cavity.

On the other hand, if substructures are absent, then dust drift depletes most of the solid material in the protoplanetary disk, which leads to a disk with a faint emission in the continuum and a cavity driven exclusively by photoevaporation, similar to the transition disks found in the mm-faint population. We note that these analogies between smooth disks and mm-faint population, and between structured disks and mm-bright population, hold independently of the dust opacity (see Sect. 4.4).

In summary, our model suggests a synergy between photoe-vaporation and substructures, where photoevaporation is responsible for creating deep cavities in the gas and dust (both micron and millimeter sized components), while the presence or absence of early substructures determines the dust trapping and the disk millimeter flux, but not necessarily the observed cavity size. Without photoevaporation a moderate sized planet would not be able to create a deep cavity, and without early substructures photoevaporation would not be able to trap the amount of dust observed in bright transition disks.

Regarding the absolute dust mass and millimeter flux values found in our simulations for the smooth and structured disk models, we highlight that these depend on the initial dust mass used in the simulation setup, and it is more meaningful to look at the fraction of dust that was retained or lost relative to the initial value. For our fiducial setup that means that ≈25% of the initial dust mass was still present for the structured disk, and only ≈2% for the smooth disk. Transition disks with higher millimeter fluxes than those found in our current paper, such as LkCa 15 or GM Aur (Facchini et al. 2020; Huang et al. 2020, respectively), could be explained if we consider, for example, a disk that is initially more massive.

Another characteristic of transition disks are their relatively high accretion rates. Previously these seemed to be incompatible with photoevaporation models, which typically overpredicted the fraction of non-accreting disks with large cavities (Owen et al. 2011; Picogna et al. 2019). However, the works of Morishima (2012) and Gárate et al. (2021) showed already that a dead zone in the inner regions leads to long-lived inner disks, capable of sustaining high accretion rates while photoevaporation opens a cavity in the outer disk, and the work of Ercolano et al. (2018) and Wölfer et al. (2019) suggests the fraction of accreting disks with photoevaporative cavities could be higher if the disks are relatively depleted in carbon and oxygen.

Based on the measurements of X-ray luminosities (Preibisch et al. 2005), on the models of dead zones that suggest that these should be both common and stable (Delage et al. 2022, and in prep.), and on the abundance of observed substructures in disks (e.g., Williams & Cieza 2011; Andrews et al. 2018), which additionally are expected to form early in the disk evolution (Stadler et al. 2022), we can expect for all of the above-mentioned ingredients to influence evolution of protoplanetary disks. Then, in a combined disk model we would have that photoevaporation takes care of opening wide cavities, dead zone properties set the lifetime of the inner disk and accretion rate onto the star, and the dust trapping in early substructures (by a Saturn-to-Jupiter mass planet for example) determines the disk flux.

For example, our work shows that bright transition disks such as DM Tau, could still be explained through photoevaporation (see Sect. 6.3 by Francis et al. 2022, though a higher initial disk mass would be most likely required to explain this particular object), provided that there are additional dust trapping mechanisms acting during the early stages of disk evolution to explain the millimeter fluxes. With this new information taken into account, we can relax the constraints on massive planets as the sole cause of transition disks with wide gaps, and also relax the requirements on planet migration to link the transition disk and exoplanet populations (van der Marel & Mulders 2021).

5.2 The degeneracy of planet properties when photoevaporation is considered

In this section we want to emphasize an (unfortunate) consequence suggested by our model in terms of the characterization of planet candidates within transition disk cavities. Our model shows a disk in which photoevaporative dispersal controls the size of the cavity that expands from the inside out, and primordial substructures (which in this subsection we assume to be caused by planets) control the dust trapping and millimeter flux. Here, we can distinguish two stages in the photoevaporative dispersal process. An early stage when the photoevaporative cavity is smaller than the planet orbit (which leads to the formation of two dust rings), and a later stage where the photoevaporative cavity size is larger than the planet orbit. In the latter scenario we find that the size of the cavity is completely independent of the planet orbital location, and the disk flux in the millimeter continuum is only mildly sensitive to the gap amplitude and location (Figs. 9 and 10).

Therefore, our model suggest that there is to a strong degeneracy in which we cannot draw reliable constraints on the mass or location of a planet candidate within a transition disk based on the properties of the cavity. A wide cavity with a bright ring in the millimeter continuum might be very well caused by a superJupiter mass planet near the cavity edge, or by photoevaporative dispersal and a Saturn mass planet hidden in the inner regions. In particular, a planet inside a photoevaporative cavity would be effectively disconnected from the rest of the protoplanetary disk, and would not display a detectable circumplanetary disk during the dispersal process, and this could in principle explain why we have not detected more planet companions or circumplane-tary disks other than PDS70 (Keppler et al. 2018; Benisty et al. 2021).

To further illustrate this situation, we show a synthetic observation of our structured disk model in Fig. 11, which was generated using the SIMIO⁶ package (Kurtovic et al., in prep.) to post-process our radiative transfer models as if they were observed by ALMA with the template configuration of Elias24 (Huang et al. 2018). The image sequence shows that once the photoevaporative cavity opens (1.3 Myr, left panel) two rings can be initially identified, where the inner one is caused by the dust trapped at the photoevaporative cavity edge, while the outer one is caused by the dust trapped outside the planet gap. In our simulations, the inner disk is fainter than the outer one (at this early stage) because the photoevaporative cavity was only able to trap the remaining material from the inner regions, while the primordial gap trapped most of the material available from the outer regions. However, this stage is short-lived and lasts only for ~ 0.2 Myr in our the fiducial setup. As time passes, the pho-toevaporative cavity expands, which causes the inner dust ring to merge with the outer one (1.6 Myr, middle panel), and then to continue to expand well beyond the planet orbital location (1.9 Myr, right panel). Figure 11 suggests that by the last snapshot it would be impossible to infer anything about the planet location (that was responsible for the dust trapping) from the ring morphology and location in the millimeter continuum alone.

In order to distinguish if a cavity could be carved by photoe-vaporation or by a massive planetary companion, other types of signatures should be considered. Recent theoretical studies have focused on modeling the observational signatures from photoe-vaporative winds, including the expected dust content entrained with the gas (Franz et al. 2022a,b; Rodenkirch & Dullemond 2022), which could be used to point toward the photoevaporative origin of some cavities in transition disks where a planet companion has not yet been found, though these models are strongly dependent on the dust reservoir at the cavity edge. As shown in Fig. 3, the peak of Σ_d is always at larger radii than that of Σg, which means the dust signature in the wind may be even fainter than predicted by Franz et al. (2022a); and further studies investigate the correlation between the gas and dust distributions at the cavity edge in more detail (Picogna et al. 2023).

Local perturbations in the gas kinematics, for example, can be linked to a planetary companion embedded in the gap (e.g., Perez et al. 2015; Pinte et al. 2019; Izquierdo et al. 2022). Asymmetries such as spirals would not be caused by photoevaporation, but other processes in addition to planetary companions can also cause them, such as self-gravitating instabilities (e.g., Lodato & Rice 2005; Meru et al. 2017) and shadows cast by an inclined inner disk (Montesinos & Cuello 2018); Cuello et al. 2019).

Characterizing the gas content inside the millimeter continuum cavities can also help to differentiate between the different scenarios (van der Marel et al. 2016a), since planets tend to carve deeper cavities in the dust than in the gas, while photoevapora-tion carves deep cavities in both of the components (see reviews by Owen 2016; Ercolano & Pascucci 2017). However, we note that the presence of a dead zone may also result in a long-lived inner disk inside the photoevaporative cavity (Morishima 2012; Gárate et al. 2021). This inner disk would be rich in gas and poor in dust, and would act as an accretion reservoir for the star while photoevaporation opens a cavity in the outer regions, effectively mimicking some of the features produced by a planet-carved gap.

With this section we have highlighted how different mechanisms can interact together to produce similar features to those of a very massive planet, and that this degeneracy should be taken into account in the cases where no further evidence of the expected planet candidates is found.

Fig. 11 Synthetic ALMA observations at 1.3 mm of the structured disk model at 1.3 Myr, 1.6 Myr, and 1.8 Myr, generated using the SIMIO package (Kurtovic et al., in prep.) to post-process our radiative transfer model. The image shows how our disk would look if it was observed with the same ALMA configuration of Elias24 (Huang et al. 2018, from the DSHARP sample), assuming a distance of 139 pc and an inclination of 29°. The beam size is plotted in the lower-left corner. The orbit of the primordial gap (rgap = 40 AU) is marked with a dashed line. We note that this is not intended to be a comparison “with” Elias24.

5.3 On the relic disk problem

Several explanations have been proposed to solve the overpre-diction of relic disks by photoevaporative models. Some studies focus on reducing the predicted fraction through faster dispersal processes, such as thermal sweeping (Owen et al. 2013), and low carbon and oxygen abundances (Ercolano et al. 2018; Wölfer et al. 2019), or through long lived inner disks with dead zones that can sustain the high accretion rates for longer times (Gárate et al. 2021). Another possibility is that the radiation pressure is very efficient at removing dust during photoevaporation dispersal, since the stellar photons can transfer their momentum directly to the solid particles, which can lead to very faint disks that would be hard to detect in the IR (Owen & Kollmeier 2019).

In this work we find that the SEDs from Fig. 7 still display a high amount of FIR emission, which is above the median of class II disks from nearby star-forming regions (Ribas et al. 2017), and in line with the high FIR luminosities found in transition disks (Espaillat et al. 2014). However, we note that these emissions are only possible in our model because the dust removal through entrainment with photoevaporative winds is very inefficient, leading to low dust-loss rates, even after considering different entrainment parameters (see Appendix A and Hutchison et al. 2016; Hutchison & Clarke 2021; Booth & Clarke 2021). If we considered the additional effect of radiation pressure and magneto-thermal wind models to remove solid material (Owen & Kollmeier 2019); Rodenkirch & Dullemond 2022, respectively), we would likely find also a deficit of FIR that is more in line with the predictions for relic disks.

Another mechanism that could reduce the disk flux is the conversion of dust particles into planetesimals by streaming instability (Youdin & Goodman 2005), which tends to limit the optical thickness at the dust ring to values of τ_v ≈ 0.5 as shown by Stammler et al. (2019), however this would occur only in the cases where the local dust-to-gas ratio is high with ϵ ≳ 1, and would stop as soon as the dust content drops below this threshold. Because of the self-regulating nature of planetesimal formation, we do not expect it to significantly reduce the millimeter flux. We further note that a recent study of Carrera & Simon (2022) questions the efficiency of the streaming instability as a planetesimal formation scenario in pressure bumps, which would also reinforce our claim.

With all these different mechanisms at play it becomes unclear whether we should expect photoevaporating disks to be bright or faint in the FIR. We propose that the answer depends on the presence or absence of an inner gas disk during the dispersal process, such as the ones sustained by dead zones (Morishima 2012; Gárate et al. 2021). In the case that an inner disk is absent, the edge of the photoevaporative cavity (where most of the dust is trapped) should be directly irradiated by the central star, leading to higher photoevaporation rates, and efficient dust clearing by radiation pressure (Picogna et al. 2019; Owen & Kollmeier 2019). Instead, if an inner disk is present, it could cast a shadow on the outer disk (see Ueda et al. 2019) and shield the edge of the photoevaporative cavity from the direct irradiation, slowing its dispersal and reducing the dust loss rates found in Owen & Kollmeier (2019).

Testing this idea however, goes beyond the scope of this work, and requires a proper consideration of the dust distribution of the inner disk in the presence of a dead zone and photoevap-oration, the effect of accretion heating which may increase the scale height of the inner disk, and resolving the inner edge of the dead zone where large amounts of dust can be trapped (Ueda et al. 2019). If the presence of an inner disk during dispersal is correlated with the disk mass (which is likely if the inner disk is sustained by a dead zone, see Turner et al. 2007; Delage et al. 2022), then we would expect that more massive disks are more likely to become accreting transition disks that retain a bright dust component, while less massive disks become relic disks that both disperse quickly and lose their dust content due to radiation pressure. We expect that only vertically thicker inner disks will be able to block enough stellar irradiation to slow down the dust removal from the outer disk, and it remains to be tested whether or not such inner disk scale heights are achievable.

Our preliminary results in Appendix A on this aspect are however inconclusive, as we find that even a perfect entrainment scenario is not enough to completely remove the FIR excess. While dust removal by radiation pressure was shown to reduce the FIR emission by Owen & Kollmeier (2019, see their Fig. 11), we find that it leads to dust loss rates that are similar to those from our current work for our fiducial parameters. Thus, our model seems unable to create the undetectable relic disks so far.

Other avenues to proceed could be to reconsider the relic disk problem in the light of recent results, such as the non-membership of some of the systems considered in the Hardy et al. (2015) study (see Galli et al. 2020; Luhman 2020; Michel et al. 2021), and that some of the faint transition disks mentioned in Owen et al. (2012) might be actually inclined disks (van der Marel et al. 2022), before drawing further conclusions regarding the nature of relic disks.

Finally, we note that two recent observations of the disks J16090141 and J16070384 by van der Marel et al. (2022) could be cataloged as relic disks candidates, since they show most of the expected features, that is, large cavities in the millimeter continuum, low millimeter fluxes, low accretion rates (Alcalá et al. 2014, 2017), and low emission in the FIR (Ansdell et al. 2018), though it is still necessary to better characterize the gas depletion of the inner disk in objects before determining if they are compatible with photoevaporative dispersal scenario, and also to consider that these disks are orbiting around low mass stars.

6 Summary

In this work we performed numerical simulations to study the effect that primordial substructures have on the dust evolution and millimeter flux of protoplanetary disks undergoing photo-evaporative dispersal. Our simulations show that the presence of a primordial substructure, specifically in terms of the dust trapping efficiency, determines the flux that the disk displays during its dispersal. Once photoevaporation opens a cavity in the inner regions, further dust loss due to drift is prevented, with dust removal due to wind entrainment having little impact on the evolution of solids. Therefore, disks that developed early substructures are bright when observed in the millimeter continuum, while disks that were smooth during its early evolution are faint, with both types of disks maintaining a relatively constant flux as the photoevaporative cavity expands (Fig. 6).

From our parameter space exploration, we learn that in order to have a bright disk, while undergoing photoevaporative dispersal, the main requirement is that an early substructure with an amplitude of A_gap ≳ 2 exists (approximately what a Saturn mass planet would cause). Substructures located closer to the star could trap more dust and lead to brighter disks, while higher X-ray luminosities that lead to an earlier dispersal also reduce the amount of dust lost due to diffusion and drift prior to the cavity opening.

The millimeter fluxes calculated for the smooth and structured disk models are respectively comparable to those measured in the mm-faint (F_mm < 30 mJy) and mm-bright (F_mm > 30 mJy) populations of transition disks, a result that holds for two different opacity prescriptions. It is possible then, that at least some observations of transition disks correspond to dispersing disks where early substructures were present (bright disks) or absent (faint disks).

While it is not new that the presence of substructures determines the flux of a disk or that photoevaporation leads to the inside-out expansion of a cavity, the combination of both results has a much more distinct consequence: Bright transition disks with a large cavity are not necessarily caused by massive superJupiter mass planets. Instead, these disks can also be created by smaller planets, with masses within the range of Saturn and Jupiter that trap the dusty material at earlier times (without strong constraints on their location), and a photoevaporative dispersal process, where the edge of the expanding cavity drags all the remaining solids in the protoplanetary disk to larger radii. This scenario could also explain why we have not detected more planets in transition disk observations: they might be both less massive, and located further inside their cavities than previously expected.

From the predicted SEDs, we find that our photoevaporating disk models display a bright FIR signal, which could be problematic from the theoretical point of view, since there are no detections of non-accreting transition disks (i.e., relic disks) with such high FIR fluxes. Studies considering dust removal by radiation pressure (Owen & Kollmeier 2019), and long-lived inner disks Gárate et al. (2021) that block a fraction of the stellar flux received by the outer disk can alleviate this discrepancy between the current model and observations. Alternatively, surveys that focus on the NIR or MIR emission to detect protoplanetary disks could have missed these transition disks that would display only an FIR component.

The simulations presented in this paper, along with the previous work from Gárate et al. (2021), suggest that a comprehensive disk evolution model could explain the observed properties of transition disks, where photoevaporative dispersal is responsible for opening wide cavities, substructures created by planets of moderate mass can trap the dust required to produce the fluxes measured in the millimeter continuum, and dead zones in the inner regions lead long-lived inner disks capable of sustaining the measured accretion rates. The versatility of a combined model opens several pathways to explain the observations of transition disks, and relaxes the constraints imposed on models that rely exclusively on giant planets to explain all the features.

To conclude, we highlight the relevance of the synergy between multiple ingredients in disk evolution, in this case photoevaporation and gap-like substructures, which when considered together can explain wider range of the observed features in protoplanetary disks, with fewer constrains.

Acknowledgements

We would like to thank the anonymous referees for their reports, insights, and suggestions that greatly improved the extent of this paper. We also thank Jochen Stadler for the help with the RADMC3D setup, Nienke van der Marel and James Owen for the insightful discussions, and for sharing a selection of the SEDs and mass loss rates profiles, respectively, that were used for comparison in this paper. The authors 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, from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement No 714769, by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) Ref no. FOR 2634/2 (ER 685/7-1, ER 685/8-2), and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311.

Appendix A Dust entrainment prescriptions

Appendix A.1 Entrainment parameters

For this work, we used a prescription of dust entrainment that removes dust grains smaller than a certain size a_w, located above a certain scale height h_w. The values used and the mass loss rates obtained in this work are consistent with the values found in the studies of dust entrainment by Franz et al. (2020) and Franz et al. (2022a). However, our implementation is still a simplification of the of dust dynamics in the wind launching region, as the maximum entrainment size and scale height of the launching wind region change with radius. To address this limitation, we perform additional simulations with varying the two parameters a_w and h_w, and compare the evolution of the dust loss rate over time. The remaining parameters are the same as in the fiducial smooth disk simulation.

Figure A.1 shows that neither the variations in the entrainment scale height, nor in the maximum entrainment size have a major impact on the dust loss rate evolution, which is expected, since most of the mass loss contribution comes from the smaller grain sizes. Only when the entrainment scale height is h_w = 1 h_g we observe a constant dust loss rate toward the end of the disk lifetime, however we note that this value is highly unlikely, since Franz et al. (2020) and Rodenkirch & Dullemond (2022) show that the wind launching region is located between 3.5h_g and 5.5h_g. Based on this parameter study, we do not expect for the exact values of the entrainment prescription to influence our results.

Fig. A.1 Evolution of the dust loss rate for the fiducial smooth disk model, using different entrainment parameters aw (top) and hw (bottom).

Appendix A.2 Comparison with Booth and Clarke (2021) entrainment

The study of Booth & Clarke (2021) shows that solids can be lifted to the disk surface due to an efficient advection triggered by the photoevaporative mass loss, which may make the dust grains more easily entrained with the photoevaporative wind.

To compare how their prescription, which defines the effective size for dust grains are entrained would affect our model, we run two additional simulations based on our fiducial setup (with and without early dust traps), implementing their equation 25: $$\begin{array}{c}{a}_{\text{crit}}=0.63\left(\frac{{{\displaystyle \dot{\text{Σ}}}}_{\text{w}}}{{10}^{-12}\text{g\hspace{0.17em}}{\text{cm}}^{-2}{\text{s}}^{-1}}\right){\left(\frac{r}{10\text{AU}}\right)}^{3/2}\\ \times {\left(\frac{z_{\text{IF}}}{4H_{\text{IF}}}\right)}^{-1}{\left(\frac{M_{\ast }}{M_{\odot }}\right)}^{-1/2}{\left(\frac{{\rho }_s}{1{\text{gcm}}^{-3}}\right)}^{-1}\text{μm,}\end{array}$$(A.1)

assuming a launching height of z_IF = 4H_IF. This prescription then replaces our dust-to-gas loss ratio expression described in Equation 13.

The evolution of the mass and mass loss rate using the Booth & Clarke (2021) prescription is shown in Figure A.2, where we do not notice almost any differences from our fiducial models in terms of total mass of solids, and observe that the mass loss rates derived from their more accurate entrainment prescription are even lower than those considered in this work, making them negligible for the global disk evolution.

Fig. A.2 Comparison between the evolution of the dust mass (top) and dust loss rate (bottom), for the fiducial disk models (black) against the entrainment prescription by Booth & Clarke (2021, gray).

The lower mass loss rates can be explained by looking at the entrainment sizes derived from the Booth & Clarke (2021) prescription, which are typically lower than those of our fiducial model, as shown in Figure A.3.

Fig. A.3 Dust entertainment size given by the Booth & Clarke (2021) prescription, for the standard photoevaporative mass loss rate (solid) and the open cavity mass loss rates (dashed), for the fiducial X-ray luminosity of Lx = 1030 erg s−1.

Appendix A.3 Maximum dust entrainment

A limit test of dust entrainment in photoevaporative winds would be the case in which the dust loss rate is always proportional to the local (vertically integrated) dust-to-gas ratio, with ${\dot{\text{Σ}}}_{\text{d}}=ϵ{\dot{\text{Σ}}}_{\text{g}}$, which ignores the effects of growth and settling. This prescription would lead to the highest possible dust loss rate from entrainment alone, and to the faintest emission from the dust in IR that can be obtained from this model.

Figure A.4 shows a comparison between the resulting SED from our fiducial smooth disk simulation when considering the standard entrainment prescription (Equation 13), against the limit case when ϵ_w = ϵ. We see that the emission peaks between 40 µm and 70 µm, with a value of F_IR ≈ 10¹⁰ erg cm⁻² s⁻¹ when the ϵ_w = ϵ. prescription is considered, which is approximately half of the emission from the standard entrainment prescription (Equation 13).

Fig. A.4 SED for the smooth disk model with two different prescriptions for the dust entrainment ϵw The snapshots are taken when the cavity size is approximately rcavity = 100 AU, assuming a distance of d = 140 pc, and no inclination.

Nevertheless, the IR emission is still higher than the median of class II objects from Ribas et al. (2017), and therefore should also be detectable with current instrumentation, despite the enhanced dust entrainment and loss rates. In other words, this object could not represent one of the yet un-detected relic disks (Owen 2016).

In order to compare our mass loss rates with those found by Owen & Kollmeier (2019), which include dust removal by radiation pressure, we map the location of the photoevaporative dust trap and dust surface density into the mass loss grid shown in Figure A.5 (top panel), and extract the corresponding mass loss rate of solids. We find the expected dust loss rates by radiation pressure are comparable or lower than those found in our fiducial model, and we do not expect them to remove a significant amount of dust for our current setup. Higher dust loss rates could be reached if the dust surface densities were lower, as the mass loss grid reaches its maximum at values of Σ_dust ~ 3 × 10⁻³ g cm⁻², however it seems unlikely that the additional contribution of radiation pressure alone can remove enough dust to completely dim the FIR emission, given our radiative transfer calculations of the SED emission shown in Figure A.4.

Fig. A.5 Top: Mass loss rates from Owen & Kollmeier (2019, Figure 6) due to radiation pressure, with the photoevaporative trap location and dust surface density evolution of our fiducial model plotted on top (dashed line). Bottom: Mass loss rates for the fiducial smooth disk model (solid), the disk model with uniform entrainment (ϵw = ϵ, dotted), and the estimated loss rates from the Owen & Kollmeier (2019) radiation pressure (dashed), given our fiducial model properties.

References

All Tables

All Figures

	Fig. 1 Example of the gas loss rate profile, following the photoevaporation model from Picogna et al. (2019, their Eqs. (2) and (4)), for Lx = 1030 erg s−1. The figure shows the mass loss rate before the photoevap-orative cavity opens (solid line) and after the photoevaporative cavity opens (dashed line), with the cavity edge located at rcavity ≈ 30 AU after 1.9 Myr of evolution.
In the text

	Fig. 2 Top: Mass evolution of the disk mass in gas (red lines) and dust (blue lines). Bottom: Evolution of the mass loss rate of gas and dust by photoevaporative winds. The markers indicate the moment when photo-evaporation opens a cavity in the inner disk (“+” for the initially smooth disk, “x” for the initially structured disk).
In the text

	Fig. 3 Evolution of the surface density profiles disks from t = 1 Myr, and plotted every 0.1 Myr (solid lines, with line opacity increasing with time). The dust surface density accounts for all the grain sizes. The initial condition is shown with dashed lines.
In the text

	Fig. 4 Dust size distribution for a smooth disk at 1.2, 1.6, and 2.1 Myr. The drift (cyan) and fragmentation (pink) growth limits are also indicated.
In the text

	Fig. 5 Dust size distribution for a structured disk at 1.1, 1.3, and 1.6 Myr. The drift (cyan) and fragmentation (pink) growth limits are also indicated. The last two snapshots were selected to match the ones of Fig. 4 in terms of the photoevaporative cavity size.
In the text

Fig. 6 Millimeter fluxes Fmn, as a function of the size of the photoevaporative cavity for the smooth (dotted) and structured (solid) disk models. The solid and dotted lines represent flux from the optically thin approximation (Eq. (15)). The markers are the fluxes obtained with RADMC-3D (“+” for the smooth disk, “x” for the structured disk). The disks are assumed to be at 140 pc, and no inclination. The cavity size measurement is based on the dust distribution of millimeter sized grains.

In the text

	Fig. 7 SEDs for the smooth (dotted) and structured (solid) disk models, when the photoevaporative cavity size is rcavity ≈ 15 AU (gray), and rcavity ≈ 100 AU (black), assuming a distance of 140 pc and face-on. The data points show the SED of SzCha from van der Marel et al. (2016b) survey, re-scaled to a distance of 140 pc.
In the text

	Fig. 8 Evolution of the dust mass (top) and disk flux at λ = 1.3 mm (bottom, assuming a distance of 140 pc), for different X-ray luminosities Lx, with the black line corresponding to the fiducial value. The markers indicate the moment when photoevaporation opens a cavity in the inner disk (“+” for the smooth disk, “x” for the structured disk).
In the text

	Fig. 9 Same as Fig. 8, but for structured disks with different gap amplitudes, with a fixed location at rgap = 40 AU and X-ray luminosity of Lx = 1030 erg s−1. Notice that the axis scales are different from Fig. 8.
In the text

	Fig. 10 Same as Fig. 8, but for structured disks with different gap locations, with a fixed amplitude of Agap = 4 and X-ray luminosity of Lx = 1030 erg s−1. Notice that the y-axis scale is different from Fig. 9.
In the text

Fig. 11 Synthetic ALMA observations at 1.3 mm of the structured disk model at 1.3 Myr, 1.6 Myr, and 1.8 Myr, generated using the SIMIO package (Kurtovic et al., in prep.) to post-process our radiative transfer model. The image shows how our disk would look if it was observed with the same ALMA configuration of Elias24 (Huang et al. 2018, from the DSHARP sample), assuming a distance of 139 pc and an inclination of 29°. The beam size is plotted in the lower-left corner. The orbit of the primordial gap (rgap = 40 AU) is marked with a dashed line. We note that this is not intended to be a comparison “with” Elias24.

In the text

	Fig. A.1 Evolution of the dust loss rate for the fiducial smooth disk model, using different entrainment parameters aw (top) and hw (bottom).
In the text

	Fig. A.2 Comparison between the evolution of the dust mass (top) and dust loss rate (bottom), for the fiducial disk models (black) against the entrainment prescription by Booth & Clarke (2021, gray).
In the text

	Fig. A.3 Dust entertainment size given by the Booth & Clarke (2021) prescription, for the standard photoevaporative mass loss rate (solid) and the open cavity mass loss rates (dashed), for the fiducial X-ray luminosity of Lx = 1030 erg s−1.
In the text

	Fig. A.4 SED for the smooth disk model with two different prescriptions for the dust entrainment ϵw The snapshots are taken when the cavity size is approximately rcavity = 100 AU, assuming a distance of d = 140 pc, and no inclination.
In the text

Fig. A.5 Top: Mass loss rates from Owen & Kollmeier (2019, Figure 6) due to radiation pressure, with the photoevaporative trap location and dust surface density evolution of our fiducial model plotted on top (dashed line). Bottom: Mass loss rates for the fiducial smooth disk model (solid), the disk model with uniform entrainment (ϵw = ϵ, dotted), and the estimated loss rates from the Owen & Kollmeier (2019) radiation pressure (dashed), given our fiducial model properties.

In the text