Large gaps and high accretion rates in photoevaporative transition disks with a dead zone

Observations of young stars hosting transition disks show that several of them have high accretion rates, despite their disks presenting extended cavities in their dust component. This represents a challenge for theoretical models, which struggle to reproduce both features. We explore if a disk evolution model, including a dead zone and disk dispersal by X-ray photoevaporation, can explain the high accretion rates and large gaps (or cavities) measured in transition disks. We implement a dead zone turbulence profile and a photoevaporative mass loss profile into numerical simulations of gas and dust. We perform a population synthesis study of the gas component, and obtain synthetic images and SED of the dust component through radiative transfer calculations. This model results in long lived inner disks and fast dispersing outer disks, that can reproduce both the accretion rates and gap sizes observed in transition disks. For a dead zone of turbulence $\alpha_{dz} = 10^{-4}$ and extent $r_{dz}$ = 10 AU, our population synthesis study shows that $63\%$ of our transition disks are accreting with $\dot{M}_g>10^{-11} M_\odot/yr$ after opening a gap. Among those accreting transition disks, half display accretion rates higher than $5\times10^{-10} M_\odot/yr$ . The dust component in these disks is distributed in two regions: in a compact inner disk inside the dead zone, and in a ring at the outer edge of the photoevaporative gap, which can be located between 20 AU and 100 AU. Our radiative transfer calculations show that the disk displays an inner disk and an outer ring in the millimeter continuum, a feature observed in some transition disks. A disk model considering X-ray photoevaporative dispersal in combination with dead zones can explain several of the observed properties in transition disks including: the high accretion rates, the large gaps, and long-lived inner disks at mm-emission.


Introduction
The nature of the observed structures in protoplanetary disks and their relation to the disk evolution has been a subject of study for over 30 years now. Transition disks, in particular, are disks that present a deficit in the near-infrared (NIR) and/or midinfrared (MIR) emission, while still displaying the characteristic far-infrared (FIR) excess of most protoplanetary disks (Strom et al. 1989;Skrutskie et al. 1990), and represent one of the key puzzle pieces to understand the disk evolution process (see reviews by Owen 2016;Ercolano & Pascucci 2017).
The lack of NIR/MIR emission in transition disks' spectral energy distribution (SED) is attributed to a gap in the dust component, or more precisely, to the lack of hot micron sized grains in the inner disk (see Espaillat et al. 2014, for a review). Through radiative transfer models it has been possible to measure the radial extent of the dust cavity for a wide sample of transition disks (van der Marel et al. 2016b) and observations at different wavelengths have also shown that some of these objects retain a compact dust component close to the star (e.g. Espaillat et al. 2010;Benisty et al. 2010;Olofsson et al. 2013;Matter et al. 2016;Kluska et al. 2018;Pinilla et al. 2019Pinilla et al. , 2021, demonstrating that for some transition disks the gap is not completely devoid of dust. One particularly curious and challenging feature of transition disks is that a large number of them show gas accretion signatures, with rates that can be as high asṀ g ∼ 10 −8 M yr −1 (e.g. Cieza et al. 2012;Alcalá et al. 2014;Manara et al. 2014Manara et al. , 2017, indicating that there is plenty of gaseous material close to the star, or that low density material accretes at supersonic speeds into the central star, driven by a combination of winds and magneto-hydrodynamic (MHD) processes (e.g. Wang & Goodman 2017). If the inner cavities are really rich in gas, then this raises the following question: How is it possible to create a prominent dust gap dust, while still retaining a long lived inner gas disk that can produce such accretion rates?
Classical theoretical models suggest that disk evolution occurs in two different stages: A viscous evolution stage, where the transport of angular momentum drives the accretion of gas onto the star (Lynden-Bell & Pringle 1974;Pringle 1981), and a dispersal stage, where the mass loss rate due to photoevaporation overcomes the accretion rate onto the star, creating a cavity and dispersing the disk from the inside out (Clarke et al. 2001;Alexander et al. 2006a,b;Alexander & Armitage 2007 Though photoevaporation models can easily predict large gaps that extend for tens of AU, these fail to explain the observed high accretion rates, since the gaseous inner disks are short-lived and quickly accreted once the gap opens (Owen et al. 2010(Owen et al. , 2011Picogna et al. 2019). As a consequence, these models tend to over-predict the fraction of non-accreting disks with gaps that extend beyond 20 AU, although the mechanisms that speed up the depletion of the outer disk, such as thermal sweeping and depletion of carbon and oxygen, can partially alleviate the discrepancy with observations (Owen et al. 2013;Ercolano et al. 2018;Wölfer et al. 2019).
On the other hand, planet-disk interactions allow for high accretion rates in the inner disk (Lubow & D'Angelo 2006), and also create a gap where the large dust particles are trapped at its outer edge (Pinilla et al. 2012a). Depending on the planet mass and disk properties (e.g. disk viscosity), the inner disk can be replenished of micron-sized particles from the outer disk. However, up to today planets have been only detected in the cavity of one transition disk (PDS 70, Keppler et al. 2018;Christiaens et al. 2019), while current observational capabilities should already have detected some of the planets inferred to explain the structures of some of them (e.g. Asensio-Torres et al. 2021). This questions the universality of planets as the potential origin of the transition disks like structures. Invoking multiple planets can explain a wider variety of structures (Pinilla et al. 2015), but it can also reduce the lifetime of the inner disk, leading to the same problems as photoevaporation models (Zhu et al. 2011).
Other models including grain growth and/or dead zones (i.e. regions with low turbulent viscosity, Gammie 1996), show that these can create transition disk structures (Dullemond & Dominik 2005;Birnstiel et al. 2012a;Flock et al. 2015) comparable to the ones observed at different wavelengths (e.g. Regály et al. 2012;Pinilla et al. 2016a).
In this work we aim to explain the accretion rates and gap sizes observed in transition disks by revisiting and expanding the model of Morishima (2012), which studied the evolution of disks with a dead zone undergoing photoevaporative dispersal. In their model, a protoplanetary disk is expected to evolve through the following steps: -The disk viscous evolution is driven by the turbulence created by the magneto rotational instability (MRI, Balbus & Hawley 1998). -In the inner regions the ionization fraction is low, creating a "dead zone" where the MRI is inefficient, and the turbulent viscosity is low (Gammie 1996). -Because of the difference in turbulent viscosity, the inner disk evolves slower than the outer disk. This causes the accretion rate at the outer "active" regions to decrease faster than in the inner "dead" regions. -Photoevaporation clears a gap outside the dead zone, once the accretion rate has dropped below the mass loss rate. Meanwhile, the inner disk is unaffected by photoevaporation. -This results in a long lived inner disk that acts as an accretion reservoir, and a large gap in the outer disk that continues to expand from the inside out due to photoevaporation (Morishima 2012; Bae et al. 2013).
The previous study of Morishima (2012) demonstrated that large gaps and high accretion rates can be recreated through the combined effects of dead zones and photoevaporation, however it is still unclear if this model can produce a fraction of accreting transition disks that is consistent with observations (see Hardy et al. 2015;Owen 2016), what is the predicted distribution of accretion rates and cavity sizes, and if the dust component can reproduce the characteristic features of transition disks, such as the deficit in the NIR/MIR emission.
Motivated by these open questions, we construct a series of population synthesis models, where we include the state-of-theart X-ray photoevaporation model from Picogna et al. (2019) and a parametric dead zone prescription, and implement them into the modular disk evolution code DustPy (Stammler and Birnstiel, in prep.). For each population we want to measure the gas accretion rate and the gap size distribution, compare it with the observed population of transition disks and determine which dead zone properties are more likely to produce high accretion rates with large gaps, without overestimating the fraction of nonaccreting disks.
To test if this model can produce a signal that is consistent with observations of transition disks, we also include a dust coagulation model, where we track the evolution of the grain size distribution during disk dispersal, and obtain synthetic SEDs and images at both NIR and millimeter wavelengths.
This paper is organized as follows. In Section 2 we present our disk evolution model. In Section 3 we describe our numerical setup of our simulations. In Section 4 we show the results of our population synthesis study, focusing on the gas evolution. Section 5 shows the evolution of a single simulation including the evolution of dust, along with synthetic images from a radiative transfer calculation. In Section 6 we discuss how our model relates to the observed population of transition disks, its caveats and follow up work. We summarize our results in Section 7.

Disk evolution model
For our model we consider an axisymmetric protoplanetary disk undergoing viscous accretion, and describe its gas surface density evolution in 1D using the following diffusion equation (Lüst 1952;Lynden-Bell & Pringle 1974;Pringle 1981): where Σ g is the gas surface density, r is the radial distance to the star (in cylindrical coordinates), ν is the kinematic viscosity anḋ Σ w is the mass loss rate due to photoevaporation. For the kinematic viscosity we adopt the Shakura & Sunyaev (1973) where α is a dimensionless parameter that controls the intensity of the disk turbulence, h g is the gas scale height, and c s = γk B T/µm p is the sound speed, with k B the Boltzmann constant, T the disk temperature, µ = 2.3 the mean molecular weight, m p the proton mass and γ = 1.4 the adiabatic index.

Photoevaporation model
Photoevaporation is believed to play a major role in the dispersal of protoplanetary disks. As the accretion rate decreases over time due to viscous evolution, disks naturally transition to a photoevaporating regime where a gap opens from inside out, once the mass loss rate surpasses the local gas accretion rate (Clarke et al. 2001;Alexander et al. 2006a,b).
For this study, we focus on the effects of X-ray photoevaporation (Ercolano et al. 2008(Ercolano et al. , 2009Owen et al. 2010), and implement the mass loss rate profiles from Picogna et al. (2019) into the sink termΣ w of Equation 1.
We refer the reader to the original Picogna et al. (2019, their Eqs. 2-5) for a detailed description of the mass loss rate, which is derived from 2D hydrodynamic models with radiative transfer calculations, and parameterized as a function of the stellar X-ray luminosity L x . Figure 1 shows the correspondingΣ w profile used for this work.
For reference, in this model an X-ray luminosity of L x = 10 30 erg s −1 corresponds to a total mass loss rate ofṀ w ≈ 10 −8 M yr −1 , and is overall higher than the previous photoevaporative model from Owen et al. (2010), making it easier to open a large gap in the gas.

Dead zone model
Turbulence driven by the magneto rotational instability (MRI) is one of the candidates to explain the angular momentum transport across protoplanetary disks (Balbus & Hawley 1998), and is triggered by the coupling between the charged particles in the gas phase and the magnetic field (Balbus & Hawley 1991).
In the inner regions of the disk, where the column density is higher and the ionization lower, the MRI is quenched or even shut off, creating a "dead zone" (Gammie 1996). In contrast with the MRI active regions, the turbulence in the dead zone is much lower, and is dominated by less efficient hydrodynamic instabilities, such as the vertical shear instability (VSI, Nelson et al. 2013;Flock et al. 2020;Manger et al. 2020).
where α a is the turbulence parameter of the MRI active region, α dz is the turbulence parameter of the dead zone, r dz is the dead zone radial extent, ∆r = r − r dz , and w = r dz /5 is the transition width from between the dead and the active regions; similar to the models used by Birnstiel et al. (2012a); Gárate et al. (2019). For reference, we show the profile in Figure 2.
This profile results in a fast viscous evolution for the outer disk (r > r dz ), and a slow viscous evolution for the dead zone in the inner disk (r < r dz ).
While this profile does not include the dependence on the surface density profile from other models (e.g. Kretke et al. 2009;Morishima 2012;Pinilla et al. 2016a), it makes it easier to study the impact of the dead zone extent on the population of accreting transition disks. We discuss this point in Section 6.3, and compare our approach with recent dead zone models.
Other mechanism angular momentum transport mechanisms, such as magnetic winds (Blandford & Payne 1982) could lead to different disk accretion profiles, and it might be worth considering their effect on the global gas surface density evolution for future studies (Suzuki et al. 2016).

Dust evolution
Since we want to include a radiative transfer calculation in this work, we also need to include a model for the dust dynamics, since these differ from the dynamics of the gas (Whipple 1972;Weidenschilling 1977), as revealed by observations of disks at different wavelengths (e.g. TW Hya, Andrews et al. 2016;van Boekel et al. 2017;Huang et al. 2018).
While small particles are indeed well coupled to the gas motion, larger particles can decouple from it, depending on their size. A useful quantity to characterize the level of coupling is the Stokes number, St = t stop Ω k , where t stop is the time necessary for the dust grains to couple to the gas motion due to their mutual drag force, and Ω k is the Keplerian angular velocity.
For spherical dust grains located at the disk midplane, we can write the Stokes number of a particle of radius a, and material density ρ s as: where λ mfp = (nσ H 2 ) −1 is the mean free path, with n the number density and σ H 2 = 2 × 10 −15 cm 2 the molecular cross section.
Because most of the mass is concentrated in large grains, and these tend to settle towards the midplane, Equation 5 is a convenient expression to describe the dust aerodynamic behavior.
It can also be useful to think of the Stokes number as a "dynamical grain size". From Equation 6 we see that small grains (with St 1) move with the viscous velocity of the gas, large boulders (St 1) do not move in the radial direction, and midsized pebbles (St ≈ 1) drift towards the pressure maximum at v d ≈ −ηv k .
In addition to advection, dust particles also diffuse according to the concentration gradient, with a diffusivity of D d = ν/(1+St 2 ) (Youdin & Lithwick 2007). The dust evolution is then described by the following advection-diffusion equation (Birnstiel et al. 2010): where Σ d is the dust surface density (of a particular dust species), andΣ w,d is the dust loss rate due to entrainment with the photoevaporative wind (Hutchison et al. 2016;Franz et al. 2020;Hutchison & Clarke 2021;Booth & Clarke 2021), which we describe below.
In this work we ignore the effects of the dust back-reaction onto the gas dynamics, since these become negligible at low dust-to-gas ratios ( ≤ 0.01) and over long ∼ Myr timescales (Gárate et al. 2020).

Dust settling and wind entrainment
Dust models predict that as particles grow in size, they also tend to settle towards the midplane (Dubrulle et al. 1995), a behavior that is confirmed by observations (see Villenave et al. 2019Villenave et al. , 2020, for a prominent example).
Assuming that the gas is in hydrostatic equilibrium with a characteristic scale height h g , the dust scale height can be approximated by (Youdin & Lithwick 2007): Then, the vertical structure of the gas and dust can be modeled with a Gaussian distribution (Fromang & Nelson 2009) as: Since small particles are well coupled to the gas and we expect the photoevaporative wind to be launched from the disk surface, we can estimate the dust loss rate due to wind entrainment as: where a w ,h w is the dust-to-gas ratio of small entrained particles (a ≤ a w = 10 µm) above several gas scale heights (z ≥ h w = 3h g ). We pick these limits based on the model of Franz et al. (2020), which indicate that only small particles that are several scale heights above the midplane can get entrained with the wind. While this approximation is simplistic, we find it effective because the dependency of the mass loss rate with the a w and h w parameters is very weak (Gárate et al., in prep.).

Dust growth
The last ingredients of our dust evolution model are growth and fragmentation which are necessary if we intend to model the simultaneous evolution of gas and dust (e.g., Brauer et al. 2008;Birnstiel et al. 2010;Drażkowska et al. 2019).
In this work, we consider that the dust phase consists of a distribution of particle species with different sizes, with their dynamics determined by their Stokes number (Equation 5), that can also evolve through sticking and fragmentation, as described in Birnstiel et al. (2010).
Without entering into the details of the coagulation model, we want to remark that the growth of a particle distribution is locally limited either by drift, when the drift timescale exceeds the growth timescale, or by fragmentation, when the collision velocity between particles exceeds the material fragmentation velocity v frag (Ormel & Cuzzi 2007;Brauer et al. 2008;Birnstiel et al. 2009).
Then, the maximum Stokes number that a particle can reach in each case is: with the maximum grain size given by St max = min(St frag , St drift ) (Birnstiel et al. 2010. Because the collision velocity between particles depends on the local turbulence parameter (Ormel & Cuzzi 2007), with: we can expect particles to grow into larger sizes in regions with low α (see Equation 12), such as the dead zone Pinilla et al. 2016a;Ueda et al. 2019).

Simulation setup
We perform our numerical simulations using the code DustPy 1 (Stammler and Birnstiel, in prep.), that can simulate the advection of gas and dust, along with the growth and fragmentation of multiple particle species, following the model of Birnstiel et al. (2010). The code also allows to load custom modules, which is how we implemented the photoevaporation and dead zone profiles described in Section 2. This study is separated into "gas only" simulations, that are fast to run and ideal for population synthesis studies (which are the main focus of this work) and a single "gas and dust" simulation, in which we solve the evolution of both components simultaneously to perform a radiative transfer calculation of a specific Matías Gárate et al.: Large gaps and high accretion rates in photoevaporative transition disks with a dead zone  [5,10,20] case of study, in order to compare with multi-wavelength observations.
In this section, we describe the parameters used to set up the protoplanetary disk, the numerical grid and the parameter space explored.

Disk setup
For our simulations we use a solar mass star, and a circumstellar disk with an initial mass of M disk = 0.1M . The initial gas surface density is defined with the Lynden-Bell & Pringle (1974) self-similar solution: This initial condition is defined by the characteristic radius r c at which the exponential drop begins, and a normalization factor Σ 0 , defined such that 2πrΣ g dr = M disk . For the gas temperature we assume that the disk is heated passively by the central star, and use the following profile: For all our simulations, the radial grid extends from 1 AU to 300 AU, with n r = 250 logarithmically spaced grid cells.
Additionally, in Appendix A we test the effect of the initial condition on our results, specifically, we present what would happen if the disk inner regions were in a quasi-steady state from the beginning of the simulation (i.e. a radially constant accretion profile, that decreases over time).

Population synthesis
In this work we construct 10 different disk populations, consisting of one "Control" population without a dead zone (i.e., α(r) = α a ), and 9 different "Dead Zone" populations, with varying radial extents r dz and turbulence values α dz . In Table 1 we list the parameter space explored for r dz and α dz . The active turbulence value of α a = 10 −3 is kept constant for all the simulations.
To construct the disk populations, we use a similar approach to the previous studies from Owen et al. (2011);Ercolano et al. (2018) and Picogna et al. (2019). Each population consists of 1000 simulations, where the X-ray luminosity is sampled from the Taurus luminosity distribution, with values between L x = 10 28 erg s −1 and L x = 10 31 erg s −1 (Preibisch et al. 2005;Güdel et al. 2007), and the disk characteristic size is sampled from a uniform distribution with values between r c = 20 AU and r c = 100 AU. To ensure that the populations are comparable between each other, we use the same sample of L x and r c across the different populations.
We track the evolution of each simulation for 10 Myr, or until the photoevaporation carves a gap that extends beyond 100 AU, saving a snapshot every 0.1 Myr 2 .
In particular we want to obtain the probability distribution of the gas accretion rateṀ g (measured at 1 AU) and the outer edge of the gap opened by photoevaporation including the presence of a dead zone, for each population model (see Section 4), and compare them to the observed distribution of accreting transition disks.

Model with dust evolution
We perform an additional simulation including dust evolution (as described in Section 2.3), and using the following parameters: L x = 10 31 erg s −1 , r c = 60 AU, α dz = 10 −4 , and r dz = 5 AU. The remaining parameters are the same as in Section 3.1.
For the size distribution we consider a logarithmic grid of n m = 141 particle species going from 5 × 10 −5 cm to 2.5 × 10 2 cm. The initial dust size distribution follows the Mathis et al. (1977) power-law, with a maximum size of a 0 = 1 µm, and we assume an initial dust-to-gas ratio of 0 = 0.01.
For the dust properties we assume that the fragmentation velocity is v frag = 10 m s −1 (Wada et al. 2011;Gundlach et al. 2011;Gundlach & Blum 2015), that the grains have a material density of ρ s = 1.6 g cm −3 and that these are compact, with the grain mass given by m = 4πρ s a 3 /3.
The dust then evolves according to the growth and advection model described in Section 2.3. Finally, we can use the dust size distribution to obtain the disk's SED and its brightness in both the millimeter continuum (1.3 mm) and scattered light (1.25 µm), using the radiative transfer code RADMC-3D 3 (Dullemond et al. 2012).
We restricted our analysis to a single "dust and gas" simulation for two reasons: First, because these simulations are computationally expensive to run. The second reason is that we are particularly interested in the scenario where a large gap opens on short timescales (≈ 1 Myr), while keeping the presence of a compact inner dust disk, as this has been seen in observations of transition disks at different wavelengths (e.g. Benisty et al. 2010;Olofsson et al. 2013 Rosotti et al. 2020). In order to ensure the fast gap opening, we use a high value for the X-ray luminosity (L x = 10 31 erg s −1 ) in this simulation.
We defer an in-depth study of the dust evolution to a follow up work, and limit ourselves to comment on the expected effect of the existing parameters based on the available results.

Results -Gas evolution
In this section we present the results of the population synthesis study with gas only simulations, showing that disks with dead zones are more likely to continue accreting during photoevaporative dispersal, than those without.
As a simplification, we will refer to a disk as "accreting" if the accretion rate at 1 AU isṀ g ≥ 10 −11 M yr −1 , and "nonaccreting" otherwise. We pick this limit in order to be consistent with previous works (e.g. Owen et al. 2011), and observational detection limits (Alcalá et al. 2017).

Proof of concept
To first illustrate how dead zones allow for a dispersing disk to sustain a high accretion rate we present the gas surface density Both disks were constructed following Section 3.1, with r c = 60 AU and L x = 10 30 erg s −1 . For the disk with a dead zone, we used an extent of r dz = 10 AU and turbulence of α dz = 10 −4 (see Figure 2). Figure 3 shows how the inner and outer parts of the disk with a dead zone evolves on different timescales. The outer disk is dispersed once the local accretion rate drops below the photoevaporation mass loss rate (Clarke et al. 2001), while the inner disk evolves on a longer viscous timescale (Morishima 2012).
In contrast, the inner region of the disk without a dead zone gets accreted onto the star immediately after the gap opens. Figure 4 shows how the disk with a dead zone maintains an accretion rate ofṀ g ≈ 10 −9 M yr −1 during its entire evolution, even after the inner disk is disconnected from the outer region. On the other hand, for the model without a dead zone the accretion rate drops quickly after the gap opens, turning it into a non-accreting disk.

Population synthesis models
The 2D histograms in Figure 5 show the probability distribution of the accretion rate and gap size for two different population models of transition disks (each generated from 1000 simula- tions and sampling every 0.1 Myr from the moment a gap opens, see Section 3.2).
The "Control" population, which consists of disks without a dead zone, shows that only a small fraction (< 1%) of the recorded transition disks are still accreting once the gap opens (withṀ g > 10 −11 M yr −1 ). Also, once the gap exceeds 20 AU in size, all transition disks have stopped accreting. The previous population studies of Owen et al. (2011) and Picogna et al. (2019) also show similar results.
The "Dead Zone" population, on the other hand, shows that a high fraction of transition disks are accreting, with 31% of them accreting at rates ofṀ g > 5.0 × 10 −10 M yr −1 and 63% at rates ofṀ g > 10 −11 M yr −1 . We also find that the fraction of accreting transition disks is higher in disks with smaller gaps. In particular, if we only look at transition disks with gaps sizes smaller than 50 AU, we find that 85% of them are accreting.
In Figure 5 we also over-plot the measured gap sizes and accretion rates of observed transition disk population (van der Marel et al. 2016b;Keane et al. 2014;Manara et al. 2014Manara et al. , 2016aManara et al. ,b, 2017Garcia Lopez et al. 2006;Pascucci et al. 2007;Najita et al. 2007;Spezzi et al. 2008;Merín et al. 2010;Donehew & Brittain 2011;Curran et al. 2011;Alcalá et al. 2014;Cieza et al. 2010Cieza et al. , 2012Romero et al. 2012;Follette et al. 2015), using the objects listed in Ercolano & Pascucci (2017,  Comparing our results with the observed transition disks, we find that the model including a dead zone can reproduce the measured accretion rates and gap sizes larger than 10 AU. However,  6. Evolution of the fraction of transition disks, relative to the whole disk population, over time for both the "Control" and "Dead Zone" populations shown in Figure 5. For simplicity, we identify a "transition disk" as a disk with a gap opened by photoevaporation. we also notice that this model does not predict disks with gaps smaller than the dead zone radial extent r dz , since the inner regions evolve too slowly for photoevaporation to carve a gap during the early stages of disk dispersal. In Section 6.1 we discuss how disks with narrow gaps could still be predicted within our framework. In Figure 6 we show the fraction of transition disks relative to the whole disk population (i.e. transition and primordial disks) at each time. In both the "Control" and "Dead Zone" population, we find that photoevaporation opens a gap in approximately 20% of the disks by 2 Myrs, and in 90% of the disks at some point of their lifetime (i.e. within 10 Myrs), indicating that while the dead zone helps to sustain the inner disk for longer timescales, it does not affect whether or not a gap opens, and this would be solely determined by the star's X-ray luminosity and the disk global properties (e.g. mass, size, and turbulence).

Effect of the dead zone properties.
To study how the properties of the dead zone affect the population of transition disks, we conduct a parameter space study for different dead zones sizes (r dz ) and turbulence parameters (α dz ), shown in Figure 7. To complement the figure we include the fraction of transition disks that present accretion rates higher than 10 −11 M yr −1 in Table 2, for the different dead zone properties. In all cases we see a higher fraction of accreting transition disks when a dead zone is considered. Even the compact dead zone population with high turbulence (α dz = 5 × 10 −4 and r dz = 5 AU) shows ≈ 20 times more accreting disks than the "Control" population without a dead zone.
From the distribution of accretion rates and gap sizes we learn that disks with low turbulence (α dz = 10 −4 ) present the higher fraction of accreting transition disks, which is expected since the inner disk will evolve on longer viscous timescales.
The effect of the dead zone size is less straightforward, with medium sized dead zones (r dz = 10 AU) producing the highest fraction of accreting transition disks for a fixed α dz (see Table 2). From our intuition we would expect larger dead zones to have longer viscous evolution timescales, with t ν ∼ r 2 dz /ν (r dz ), and therefore to survive for longer, which is indeed what happens when we increase the dead zone size from 5 AU to 10 AU. However, when further increasing the dead zone size from 10 AU to 20 AU we find that the fraction of accreting disks decreases from 63% to 55% (for α dz = 10 −4 ). We can explain this behavior if we consider that before photoevaporation opens a gap, dead zones evolve together with the whole disk, receiving material from the outer regions which is then redistributed across the inner regions within the local viscous timescale. In mid-size dead zones the material is redistributed more efficiently than in more extended ones, resulting in a higher and spatially uniform accretion rate across the inner disk, and a longer lifetime. In contrast, for more extended dead zones the incoming material from the outer disk Additionally, we also find that the fraction of accreting transition disks can be higher, if the initial mass distribution follows a quasi-steady state solution (see Appendix A).

Results -Dust evolution
In this section, we describe the evolution of the dust component for a disk with a dead zone (following the setup from Section 3.3), and show the expected emission using radiative transfer calculations. Figure 8 shows how the dust surface density (integrated for all the grain sizes) evolves along with the gas during the disk dispersal in the case of a strong X-ray luminosity (L x = 10 31 erg s −1 ).

Dust distribution
Before the gap opening, the dust drifts from the outer disk towards the inner dead zone regions, as it looses angular momentum due to the drag force from the gas (Whipple 1972;Weidenschilling 1977).
Because of the low turbulence in the dead zone region, dust particles can grow to larger sizes and therefore drift faster, de-10 0 10 1 10 2 r (AU) 10 5 10 3 10 1 10 1 10 3 (g/cm 2 ) Gas Dust Fig. 8. Surface density evolution for gas and dust (integrated for all particle sizes) in a disk with a dead zone (as described in Section 3.3). The dashed line shows the initial condition, the solid lines show the disk evolution every 0.2 Myr. The star has an X-ray luminosity of L x = 10 31 erg s −1 , which leads to a gap opening at t ≈ 1 Myr.
pleting the inner disk of solids on shorter timescales (Equation 6 and 11).
The dust size distribution (Figure 9, top panel) shows that in the dead zone the dust particles can grow up to sizes of a ∼ Additionally, we notice that our dead zone model results in a disk where particle growth is completely limited by drift, whereas disks without dead zones would be typically dominated by drift only in the outer regions (r 10 AU for α = 10 −3 , after 1 Myr), and by fragmentation in the inner ones (Birnstiel et al. 2010. Once the photoevaporation opens a gap (approximately at 1 Myr), the inner and outer disk become disconnected in terms of mass transfer.
The inner disk continues to lose material as the dust drifts inward, though we also observe that the drift timescales increase over time. This happens because when the dust-to-gas ratio decreases, the maximum Stokes number decreases as well (Equa-tion 11), which leads to slower drift velocities (Equation 6). This behavior is in line with the arguments described in Birnstiel et al. (2012b); Powell et al. (2017), in which the lifetime of disks with grain growth limited by drift is approximately equal to the growth and drift timescales. Since the growth timescale is approximately t growth = ( Ω k ) −1 , we find that the time elapsed after the gap opening (i.e. the lifetime of the inner disk as an isolated system) determines the dust-to-gas ratio present in the inner regions.
From the dust distribution of the inner region (Figure 9, middle and bottom panels), we observe that it is dominated mostly by large cm-sized grains, and that there is a lack of micron sized dust. This happens because the collision velocities in the dead zone are low (Equation 13), and prevent the small particles to be replenished by fragmentation.
Meanwhile, the solids in the outer disk drift towards the edge of the photoevaporative gap, which is the new local pressure maximum. The dust grains form a ring, in which growth is limited by fragmentation, reaching sizes of a ≈ 1 mm. As the gap expands over time, the dust ring moves along with it (Alexander & Armitage 2007; Owen & Kollmeier 2019).

Dust emission
Using the dust size distributions obtained with DustPy in Section 5.1, we perform radiative transfer calculations using the code RADMC-3D.
We calculate the opacities for each grain size using the program OpTool 4 assuming compact spheres and a dust composition consisting of water ice (Warren & Brandt 2008), astronomical silicates (Draine 2006), troilite and refractory organics (Henning & Stognienko 1996) as stated in Birnstiel et al. (2018). To account for the problems caused by the strong forward scattering peak in the phase function of large grains and the limited spatial resolution or our model, we set the scattering opacity to an upper limit within the first 5 • of the phase function. We assumed the solar type star to be a blackbody point source and used 10 7 photons for our radiative transfer calculations. First, the 3D dust temperature structure is computed which is then used to calculate the synthetic images and SEDs with the full treatment of anisotropic scattering with polarization (Kataoka et al. 2015), as shown in Figure 10 and Figure 11, respectively. The disk is assumed to have an inclination of 20°and is observed at a distance of 140 pc, which is the mean value of the nearest star forming regions. Finally, both images were convolved using a Gaussian point-spreadfunction with full-width at half-maximum (FWHM) of 0.04" × 0.04" to mimic the current angular resolution obtained by instruments like the "Spectro-Polarimetric High-contrast Exoplanet Research"(SPHERE) at the Very Large Telescope (VLT), and the Atama Large Millimeter/sub-millimeter Array (ALMA). Figure 10 (top panel) shows the 1.3 mm continuum emission of our model after the gap opening (at 1.6 Myr). We observe a signal coming from the inner disk (r ≈ 10 AU) where the dead zone is located, and from a ring of dust located at r ≈ 40 AU, which corresponds to the pressure trap at the outer edge of the photoevaporative gap. The region in between is devoid of any emission, as expected from the dust distribution seen in Figure 9.
The total emission in the 1.3 mm continuum disk is F mm = 14.2 mJy. In terms of the classification of transition disks according to their fluxes (Owen et al. 2012), our model would be classified as a mm-faint disk (F mm ≤ 30 mJy). This is at odds with the observed properties of disks with large gaps (r gap 20 AU) and high accretion rates (Ṁ g ∼ 10 −9 M yr −1 ), which are more likely to have mm fluxes above 30 mJy (see Owen 2016, Figure  5). We discuss the implications of the low millimeter flux found in our model in Section 6.1.
The scattered light image (Figure 10, bottom panel) shows a wide ring, that extends approximately between 25 and 50 AU. There is a signal coming from the inner disk in the scattered light, which is mainly occluded by the chronograph that we assume (0.1 in size, which is typical for observations with SPHERE). Figure 11 shows the evolution of the SED. The SED presents a decrease in the IR emission between 1.0 Myr and 1.6 Myr, approximately between 5 µm and 30 µm, which coincides with the gap opening and the depletion of micron-sized particles in the inner disk (see Figure 9). This feature is characteristic of transition disks, and this object would be classified as such from its SED (Strom et al. 1989;Espaillat et al. 2014).

Comparison with the observed transition disk population
The discrepancy between the predictions of theoretical models and the observed population of transition disks has been an open problem for the last decade. Population synthesis models including photoevaporation predict that 88% to 97% of transition disks should be non-accreting disks with large gaps (Owen et al. 2011;Picogna et al. 2019), also called "relic disks". Observational constrains on the other hand, say that the fraction of these nonaccreting relics should be only around 3% (Hardy et al. 2015). In the literature, this discrepancy can also be found as the "relic disk problem".
Recent works from Ercolano & Clarke (2010); Ercolano et al. (2018); Wölfer et al. (2019) suggest that photoevaporation is more effective in disks that are depleted from carbon and oxygen, which results in higher mass loss rates and a faster dispersal of the outer disk. The outcome of this approach is a lower fraction of relics (up to 45%), which is much lower than standard photoevaporative models, and better suited to explain observations.
Meanwhile, models including dead zones and grain growth have also been able to reproduce some of the features observed in transition disks, including the decrease in the NIR emission due to inefficient fragmentation, and ring like structures Pinilla et al. 2016a), however these models also fail to clear extended gaps in the millimeter continuum, because the drift timescales of the large grains become too long, or because the particle trapping occurs beyond the dead zone outer edge.
In this work we present a different approach, which consist on combining the photoevaporative dispersal by X-ray radiation, with a dead zone in the inner disk. This way, both processes complement each other by covering their respective weak spots: the dead zone takes care of extending the lifetime of the inner disk, keeping the accretion rates high, while photoevaporation carves extended gaps in the gas and dust components (see also Morishima 2012).
The strongest point of this model is that it produces a high fraction disks with large gaps (r gap 20 AU) and high accretion rates (Ṁ g ∼ 10 −9 M yr −1 ) at the same time, predicting that up to 63% of all transition disks should be accreting. This fraction is higher than previous population synthesis estimates, though it is still not high enough to match the observational constrains by it-self (Hardy et al. 2015), since we still find that 37% of transition disks should be relics (i.e. non-accreting).
The predicted fraction of relic disks could be even lower if the effects of carbon depletion are considered, with the respective increase in the mass loss rate and faster dispersal of the outer disk (Ercolano et al. 2018), or if the disks inner regions are able to survive for longer times, by reaching a quasi-steady state during earlier stages of disk evolution for example (see Appendix A).
Alternatively, it could also be that many disks become too faint to be observed if dust is removed efficiently, either by drift during early stages of disk evolution, or by a strong coupling with the photoevaporative wind (Owen & Kollmeier 2019). While in our simulations the dust loss rate is always too low to affect the total dust mass, a more detailed modelling of the gas drag at the edge of the photoevaporative gap could yield higher dust loss rates at lower scale heights, after the dispersal of the inner disk.
In the context of the whole disk population (i.e. counting primordial and transition disks), our results seem to match the fraction of transition disks shown in Currie & Sicilia-Aguilar (2011), especially around 2 Myrs, with the dead zone having little to no influence on whether a gap is opened or not. At later times, we do predict more transition disks than those observed, but this discrepancy could also be bridged if we consider that some of our disk should become fainter over time, again due to dust evolution, entrainment and drifting.
We should keep in mind that the accuracy of our predictions depends, of course, on the validity of our dead zone model (Equation 4), since the fraction of accreting transition disks varies greatly across different dead zone parameters (see Table 2). We discuss this point in Section 6.3, where we calculate what a self-consistent dead zone model would look like.
One apparent limitation of our model is that it is unable to produce enough transition disks with narrow gaps (r gap 10 AU, see Figure 7), since by construction the combined effects of dead zones and photoevaporation favor opening a gap beyond the dead zone outer edge (Morishima 2012). This seems at odds with the observations, which indicate there are approximately as many transition disks with gaps narrower than 10 AU, as disks with wider gaps, with sizes between 10 AU and 100 AU (van der Marel et al. 2016b). However, most of the gaps in transition disks with sizes smaller than 10 AU are still unresolved by observations, and they are characterized by SED fitting. This means that the current gap sizes reflect the deficit of small dust grains in the inner regions, but may not correspond to an actual gap in the gas or in the millimeter continuum. We could still predict this type of narrow gaps within our framework, if we consider that dust growth in dead zones can deplete the inner regions of small grains, and create SED signatures that are alike to that of a transition disk ). Therefore, a disk could display a narrow gap signature in their SED due to a dead zone at earlier times, and a wide gap signature in both their SED and millimeter continuum at later stages due to the combined dead zone and photoevaporation effects described in Morishima (2012) and this work.
Besides the accretion rate and gap sizes, the observed transition disks seem to be divided in two groups according to their millimeter flux F mm (Owen et al. 2012). Disks with smaller gaps tend to be "mm-faint" (F mm < 30 mJy), while disks with larger cavities are "mm-bright" (F mm > 30 mJy), with the fluxes rescaled to a distance of 140 pc (see Owen 2016, Figures 3 and 5). Instead, the millimeter flux of our dust evolution model is only F mm = 14.2 mJy, despite presenting a large cavity and a high ac-cretion rate. Simulations with different parameters that were not included in this paper also show similar fluxes, around 10 mJy. This suggest that the dust content in our disks is relatively low, in comparison with with the observed transition disks, or that we are underestimating the dust opacity.
We can explain why our disks appear to be dust depleted by looking at our initial gas density profile (Equation 14), which is smooth and does not have any substructure. In such a disk the dust particles can drift inwards and into the star very quickly, since there are no pressure maxima to act as local dust traps (Pinilla et al. 2012b). In contrast, observations of protoplanetary disks show that these are rich in substructures such as rings and spirals (e.g. ALMA Partnership et al. 2015;Andrews et al. 2018), indicating that the underlying gas profile is not smooth. Recent models also indicate that dust traps are necessary in order to explain the size luminosity relationship observed in protoplanetary disks (Tripathi et al. 2017, Zormpas et al. in prep.). Including pressure bumps in our model at early times, such as a gap created by a planet, would help our disk to retain a higher fraction of dust, increasing its millimeter emission and reducing the discrepancy with observations (Pinilla et al. 2012a(Pinilla et al. , 2020. Despite this limitation in our model, it seems that the presence of a dead zone in the inner regions is a promising ingredient to explain the observed transition disk population. For a follow up work, we intend to quantify the millimeter continuum emission of disks undergoing photovaporative dispersal, and test if it is possible to retain enough solids to match the observed fluxes from transition disks through trapping (Pinilla et al. 2020).
Additional improvements for the model would be to include stars of different mass, considering their respective photoevaporation mass loss rates (Komaki et al. 2021), the evolution of the stellar luminosity during the first Myr after formation (Kunitomo et al. 2021), and the effect of carbon depletion, which leads to stronger photoevaporation rates (Ercolano et al. 2018;Wölfer et al. 2019).

Long-lived inner disk and comparison with observations
The large diversity in the SED morphology of transition disks (e.g., Cieza et al. 2007) provides evidence that the cavity of some of these disks is not completely depleted in dust. Some of the SEDs of transition disks exhibit a strong excess in the NIR wavelength range (e.g., Espaillat et al. 2010;van der Marel et al. 2016b), implying a compact optically thick inner disk very close to the star.
Observational efforts have been made to directly observe and characterize the inner disk of transition disks. For example, interferometric observations at the NIR have spatially resolved a compact inner disk in some transition disks (e.g. Benisty et al. 2010;Olofsson et al. 2013;Gravity Collaboration et al. 2019. Recent observations with ALMA also detect and in few cases resolve the inner disk of some transition disks, meaning that these inner disks are probably rich in millimeter-sized particles, or that there is an equivalent amount of smaller grains that lead to a similar signal in the millimeter continuum (e.g. Hendler et al. 2018;Pinilla et al. 2019;Keppler et al. 2019;Francis & van der Marel 2020;Rosotti et al. 2020).
The physical mechanism responsible for the observed dust in the inner disk of transition disks is an open question, and it is unclear how it can survive on timescales of millions of years. For the formation of the large gaps in these objects, mainly three possibilities are debated in the literature: embedded planets, photoevaporation, and dead zones. However, each of these scenarios have their own challenges to explain a long-lived inner disk as observed at different wavelengths.
On one hand, the planet-disk interaction models may allow to have partial dust filtration at the outer edge of the gap (Pinilla et al. 2012a;Weber et al. 2018), such that the inner disk can be replenished with small dust (micron-sized) from the outer disk. However, in this scenario these small particles grow very efficiently in the inner disk once they cross the gap and thus they are lost towards the star due to efficient drift. This problem can be mitigated by including the effect of the snow line, that can change the velocities for which particles fragment (Wada et al. 2011;Gundlach et al. 2011;Gundlach & Blum 2015), and hence slow down their radial drift (Birnstiel et al. 2010;Pinilla et al. 2016b;Drażkowska & Alibert 2017). These models of planets with the snow line are capable of reproducing the NIR excess seen in the SEDs of some transition disks, but they fail in reproducing a detectable inner disk at millimeter wavelengths (see Figure A.1 in Pinilla et al. 2016b).
On the other hand, models of dead zones acting alone can reproduce an inner dusty disk, but it is expected to be much fainter than the ring-like emission in the outer disk at millimeter wavelengths (see Fig. 7 and Fig. 8 in Pinilla et al. 2016a). Additionally, in these models there is not a clear depletion of the gas surface density inside the mm-cavity, which seems to be highly depleted from observations of 12 CO and its isotopologues (van der Marel et al. 2016a). To solve this problem, it was suggested that the inclusion of an MHD wind in the simulations could reduce the gas surface density in the inner regions, but would worsen the problem of keeping a long-lived inner disk.
Finally, photoevaporation models previous to this work predicted highly depleted cavities, in both gas and dust (Alexander & Armitage 2007;Owen & Kollmeier 2019), but these ignored the combined effect of a potential dead zone in the photoevaporative dispersal process.
The models presented in this paper can explain a dusty inner disk that remains millions of years, and that is bright at millimeter emission.
The maximum grain size in the inner disk is limited by drift, and therefore by the local dust-to-gas ratio, which is gradually decreasing with time. As a consequence the radial drift velocity of the particles in the inner disk also decreases, prolonging their lifetime ( Figure 9). This solves the conundrum of the origin and survival of these inner disks in transition disks.
From our synthetic observations at millimeter wavelength ( Figure 10), we find that the surface brightness in the inner disk is slightly fainter (by a factor of a few) than the ring at the outer edge of the photoevaporative gap. A similar feature has been observed in some transition disks (e.g., TCha and CIDA1, Hendler et al. 2018;Pinilla et al. 2021, respectively), where the inner disk has a comparable brightness to the outer disk. The millimeter sized particles in the inner regions extend until the outer edge of the inner disk, as given by the gas distribution, indicating that the dust can be retained for long timescales (along with the gas component) without the need of additional dust traps.

Validity of the dead zone model
In this work, we have modelled the presence of a static dead zone by using the α(r) profile described in Equation 4, which is defined by the free parameters α dz , α a , and r dz .
Our parameter space study in Section 4 shows that the fraction of accreting transition disks is heavily influenced by the turbulence in the dead zone α dz , and the dead zone extent r dz . For this reason we are interested to see if a more complex model could produce similar profiles in terms of the dead zone turbulence and radial extent. We perform a quick validation test using the framework developed in Delage et al. (2021). The authors have built a 1+1D magnetically-driven disk accretion model for the outer disk regions (r 1 AU), where the local mass and angular momentum transport are assumed to be controlled solely by the MRI and hydrodynamic instabilities such as e.g. the VSI. From such a model, we can derive an effective turbulent parameterᾱ that is self-consistently determined from detailed considerations of the MRI and the non-ideal MHD effects (Ohmic resistivity and Ambipolar diffusion) -given stellar properties (mass and luminosity), disk mass and dust properties. Their model accounts for the main physical processes happening in the outer regions of protoplanetary disks: (1) irradiation from the forming star; (2) dust settling; (3) ionization from stellar X-rays, galactic cosmic rays and the decay of short/long-lived radionuclides; (4) disk chemistry. Using their framework, we investigate two scenarios. First, the effective turbulenceᾱ is calculated from a gas surface density following the standard self-similar solution from Lynden-Bell & Pringle (1974). Since this scenario leads to a disk structure out of equilibrium (the accretion rate is radially non-constant, implying that some regions of the disk can accrete faster than others), we consider a second scenario where the effective turbulenceᾱ is derived by assuming steady-state accretion. In the latter, the gas surface density is self-consistently derived through an iterative process, alongside the effective turbulenceᾱ, to ensure the accretion rate to be radially constant. For each scenario, we consider two X-ray luminosities (10 30 and 10 31 erg s −1 ). The stellar and disk mass are the same as described in Section 3.1. We further assume a constant dust-to-gas ratio of = 0.01 as well as grains of size a = 1 µm, consistent with the setup described in Section 3.3. Figure 12 shows the α(r) profiles derived for the four simulations, using the framework described in Delage et al. (2021). Depending on the X-ray luminosity used or the scenario considered, the salient results are as follows: the dead zone extent r dz ranges from 22 AU to 40 AU, where the turbulence quickly changes from a low to high regime; the mean turbulence in the dead zone α dz ranges from 1.3 × 10 −4 to 3.7 × 10 −4 ; and the mean turbulence in the MRI-layer α a ranges from 1.8 × 10 −3 to 5.7 × 10 −3 .
Based on these calculations, we find that the dead zone radial extent should be larger than the one used in our models, which in principle would mean that the fraction of accreting transition disks could be better represented by populations with r dz = 20 AU. We also find that the mean values for the turbulence in the dead zone are covered by our parameter space exploration. Furthermore, the highest value found (α dz = 3.7 × 10 −4 ) would suggest that our choice of α dz = 5 × 10 −4 is quite marginal, and for moderate X-ray luminosities (L x = 10 30 erg s −1 ) the dead zone turbulence remains between α dz = 10 −4 and 2 × 10 −4 .
We expect that a population synthesis model using the dead zone model from (Delage et al. 2021) would also present gaps that extend beyond 20 -40 AU, and a fraction of accreting disks above 12% at least, and possibly closer to 55% (see the models with α dz = 1 × 10 −4 and 3 × 10 −4 , for r dz = 20 AU, Table 2).
Finally, we note that our choice for α a underestimates the mean values found by the more complex model. A higher turbulence in the MRI-active layer could result in the photoevaporative gap to open earlier, which we expect to result in a higher dust content in the outer disk, without decreasing the fraction of accreting disks found in this work. Indeed, α dz and r dz are the key parameters here, since these determine the lifetime and viscous timescale of the inner disk.

Summary
In this work we present a disk evolution model that includes Xray photoevaporative dispersal and a dead zone prescription, in order to explain the accretion rates and extended gap sizes observed in transition disks.
Our population synthesis study shows that considering a dead zone in the inner regions can easily result in transition disks with extended gaps (r gap 20 AU), and long lived inner disks capable of sustaining high accretion rates (Ṁ g ∼ 10 −9 M yr −1 ) during the dispersal process. This is a result of the differential evolution of the inner and outer disk. Because the outer turbulent disk evolves on shorter viscous timescales, its accretion rate decreases faster and enters into the photoevaporative dispersal regime earlier. Meanwhile, the viscous timescale of the inner disk is longer due to the low turbulence in the dead zone, allowing it to sustain a relatively high accretion rate.
For dead zone properties of α dz = 10 −4 and r dz = 10 AU, we predict that 63% of transition disks should still be accreting, withṀ g 10 −11 M yr −1 . We also find that half of these accreting transition disks display high accretion rates ofṀ g 5.0 × 10 −10 M yr −1 . This means that photoevaporative disk dispersal could, in fact, explain the observed accreting transition disks with large gaps, though it is still necessary to explain why we do not detect the predicted fraction of non-accreting disks.
While our model does not explicitly predict transition disks with narrow gaps (r dz 10 AU), we suggest that these could still be explained through other processes, such as dust growth within a dead zone during the earlier stages of disk evolution, which would then evolve into wide photoevaporative gaps at later stages.
From our dust evolution simulations, we learn that the inner disk retains a dust component consisting of large millimeter to centimeter sized grains, while the dust in the outer disk forms a ring at the edge of the photoevaporative gap, with grains sizes between a micrometer and a few millimeters. Radiative transfer calculations based on our dust distribution model show an inner and outer disk in the 1.3 mm continuum, with the inner disk being fainter by a factor of a few, and that the SED shows an emission deficit at near and mid-infrarred wavelengths, which is a characteristic feature of transition disk spectra.
While the millimeter flux in our model appears to be fainter than the ones observed in transition disks with large cavities, this could be solved by including dust traps in the outer disk (such as the ones caused by planets), or more vigorous photoevaporation rates (as those predicted for carbon depleted disks), in order to retain a higher fraction of solids during the early stages of disk evolution.
In future studies we plan to focus on the observability of photoevaporating disks by extending our analysis of the dust component, to determine if the predicted non-accreting disks could be too faint to be observed, and if additional substructures could reproduce the observed population of millimeter bright transition disks with high accretion rates.
In summary, our model including dead zones and photoevaporation is a good candidate to explain several features observed in transition disks, from the high accretion rates, to the large gaps, and even the inner disks observed in the millimeter continuum.

Appendix A: Models in quasi-steady state
In this appendix we show an additional set of population synthesis models with a different initial condition, to complement the results of section Section 4. For this work, we studied disks that are initialized with a Lynden-Bell & Pringle (1974) profile, which consist on a power law with an exponential decay. When the turbulence profile is constant, these disks evolve in a so called "quasi-steady state", which is characterized by a constant accretion rate in the inner regions, that decays uniformly over the viscous timescale. However, if a disk has a non-uniform turbulence profile, such as the one created by a dead zone (see Figure 2), this results in a variable accretion profile, that is not in local equilibrium (see the initial condition of Figure 4, top panel).
In order to test the effect of the initial condition, we repeat our population synthesis study as described in Section 3.2, but now assuming the material is initially distributed such that the inner regions evolve in quasi-steady state, using a scaled version of the Lynden-Bell & Pringle (1974) profile: The scaling factor α a /α(r) results in a profile that has a higher density in the dead zone region, and a lower density in the outer regions, when compared to the standard Lynden-Bell & Pringle (1974) solution. The normalization factor Σ 0 is again defined such that the integrated mass is M disk = 1 M . Figure A.1 shows the corresponding distribution of the gas accretion rate and gap sizes of transition disks, and Table A.1 shows the fraction accreting transition disks using the initial condition from Equation A.1. From these new results we immediately notice an increase in the abundance of transition disks accreting at high accretion rates. In particular, we find that for the model with α dz = 10 −4 and r dz = 10 AU, the predicted fraction of accreting transition disks jumped from 62.9% to 99.7%, with approximately half of those disks displaying accretion rates above 10 −9 M yr −1 . Another difference is that now, using the scaled initial condition, the populations with more extended dead zones (r dz = 20 AU) tend to display the highest fraction of accreting transition disks (for α dz ≤ 3 × 10 −4 ).
Both the overall increment in the fraction of accreting transition disks, and its systematic increase with the dead zone size r dz , can be explained if we think that the new initial condition represents disk in which the inner regions are "saturated" of material, that is, they have the required amount of gas to be in local steady state. In this scenario, the inner disk has the maximum possible mass by the time photoevaporation opens a gap, and therefore it can survive for the longest possible time.
Of course, the question now would be if the quasi-steady state initial condition described by Equation A.1 can be met. From our results in Figure 3, we see that in the presence of a dead zone, a disk initialized with a Lynden-Bell & Pringle (1974) profile slowly evolves into a quasi-steady state profile, like the one described by Equation A.1 (though the exact shape might vary), but with a lower normalization factor Σ 0 and larger characteristic radius r c , since some material must have been accreted, and the outer regions have expanded due to the transport of angular momentum.
We expect that disks with short viscous timescales in their inner regions would reach the quasi-steady state solution more easily. Accordingly, we find that the disk populations with smaller dead zones and higher turbulence parameters are be less affected by the choice of initial condition (see Table 2 and Table A.1), since they have enough time to reach the quasi-steady state regime before the photoevaporation opens a gap and disconnects the inner and outer disk from each other. In contrast, disks with extended dead zones are more affected by the choice of initial condition, since their viscous timescales are longer (of several million years).
Other factor that should affect whether a disk reaches the quasi-state state solution is the disk mass, since a more massive disk would be able to transport more material into the inner regions before photoevaporation becomes relevant for its global evolution, though the corresponding normalization factor would be lower, depending on how much material was accreted into the star.
Finally, one reason why it might be preferable to draw our conclusions from the standard Lynden-Bell & Pringle (1974) profile, instead from the quasi-steady state profile in Equation A.1, is that the dead zones could be subject to thermal and gravitational instabilities during the early stages of disk evolution (Martin & Lubow 2011), where the turbulence is cyclically reactivated in the inner regions, resulting in accretion outbursts, and for the accumulated material to quickly flushed into the star (Kretke et al. 2009;Audard et al. 2014;Gárate et al. 2019). With this point in mind, we can think of the results based on the quasisteady state initial condition as an upper limit to the fraction of accreting transition disks.
To conclude, we also show the effect of the quasi-steady state initial condition on the dust distribution and millimeter continuum synthetic observations in Figure A.2 and Figure A.3, also using an X-ray luminosity of L x = 10 31 ergs −1 , α dz = 10 −4 , and r dz = 10 AU, but taking the snapshot at 1.2 Myr, since the gap opens earlier when using the initial condition from Equation A.1. The main difference is that now the surface brightness of the inner disk is comparable to that of the outer ring in the millimeter continuum, as observed in some transition disks (Hendler et al. 2018;Pinilla et al. 2021