© The Authors 2026

1 Introduction

Planet formation remains one of the most challenging problems in modern astrophysics (Safronov 1972; Lissauer 1993; Pollack et al. 1996). Over the past three decades, exoplanet surveys using radial velocity, transit, microlensing, and direct imaging techniques have demonstrated that planetary systems are a common outcome of star formation, with the population of detected planets spanning a wide range of masses and orbital architectures (Cassan et al. 2012; Batalha et al. 2013; Mulders et al. 2018; Zhu & Dong 2021). The diversity of these systems – from compact super-Earth chains to distant gas giants – provides crucial constraints on the physics of protoplanetary disk evolution and planet formation (see reviews of Liu & Ji 2020; Drążkowska et al. 2023).

The core accretion scenario offers a broadly successful framework for linking these observations to the final planetary architectures (Pollack et al. 1996; Ida & Lin 2004; Mordasini et al. 2012; Bitsch et al. 2015; Liu et al. 2019; Liu & Ji 2020; Emsenhuber et al. 2021; Pan et al. 2025). In this picture, micron-sized dust particles coagulate into millimeter-centimeter pebbles (Blum & Wurm 2008; Birnstiel et al. 2012; Pinilla et al. 2012; Birnstiel et al. 2016), but further growth is uncertain: laboratory experiments (Blum & Münch 1993; Blum & Wurm 2000, 2008; Güttler et al. 2010) and numerical simulations (Brauer et al. 2008; Zsom et al. 2010; Okuzumi et al. 2012; Kataoka et al. 2013; Drążkowska et al. 2016) all indicate that bouncing, fragmentation, and rapid inward drift prevent pebbles from directly forming kilometer-sized planetesimals. These barriers are most severe inside ice lines or in disks with low dust-to-gas ratios, where pebbles drift inward before significant growth can occur (Drążkowska et al. 2023). The conversion of millimeter-sized solids into planetesimals therefore remains one of the least understood stages of planet formation (Johansen et al. 2014; Liu & Ji 2020; Drążkowska et al. 2023), demanding an additional, efficient concentration mechanism (Birnstiel et al. 2016; Drążkowska et al. 2016).

The streaming instability (SI) is widely considered a leading candidate for overcoming these barriers (Youdin & Goodman 2005; Johansen et al. 2007; Bai & Stone 2010; Carrera et al. 2015; Li & Youdin 2021; Lim et al. 2024, 2025a,b). By tapping the relative motion between gas and solids, the SI concentrates drifting solids into dense filaments that can collapse gravitationally into planetesimals. Its activation, however, is sensitive to the local dust-to-gas ratio, particles’ Stokes number, and disk turbulence (Bai & Stone 2010; Yang et al. 2017; Drążkowska et al. 2014; Lim et al. 2024). A local dust-to-gas ratio of a few times 10⁻² is typically required (Johansen et al. 2009b; Carrera et al. 2015; Drążkowska et al. 2016; Liu & Ji 2020). This motivates the question of how and where protoplanetary disks can enhance the dust-to-gas ratio sufficiently to trigger the SI.

High-resolution imaging of protoplanetary disks has revealed ubiquitous substructures – such as rings, gaps, and cavities – which are interpreted as pieces of evidence that dust growth and radial drift are already actively reshaping disks at very young ages (Huang et al. 2018; Dullemond et al. 2018; Andrews 2020). Pressure bumps are thought to be natural sites for halting dust drift and locally boosting the metallicity. They may arise from gaps opened by massive planets (Kobayashi et al. 2012; Zhu et al. 2012; Gonzalez et al. 2015; Eriksson et al. 2020), at the inner edges of magnetically inactive dead zones (Varnière & Tagger 2006; Lyra et al. 2009; Dzyurkevich et al. 2010), from zonal flows generated by hydrodynamical or magnetic instabilities (Johansen et al. 2009a; Dittrich et al. 2013; Bai 2014; Raettig et al. 2015; Xu & Bai 2022; Lehmann & Lin 2025), in the ice line regions (Ida & Guillot 2016; Drążkowska & Alibert 2017; Schoonenberg & Ormel 2017), or at the outer edges of photoevaporative cavities (Carrera et al. 2017; Ercolano et al. 2017). Such photoevaporative cavity edges may also influence planet migration through the associated positive density gradients, leading to outward “rebound” migration (Liu et al. 2025).

High-energy radiation, particularly X-rays from young stellar hosts, plays a critical role in clearing gas from protoplanetary disks during their late evolutionary stages (Owen et al. 2010, 2012; Ercolano et al. 2017; Wang & Goodman 2017; Jennings et al. 2018; Picogna et al. 2019). This radiation heats and ionizes the upper layers of disks, driving the photoevaporative winds that steadily deplete the gas. This process initially opens a cavity at orbital radii of a few au for a solar-mass star (closer in for a lower-mass star). As depletion continues, the cavity expands outward (Owen et al. 2010; Alexander et al. 2014; Ercolano & Pascucci 2017; Picogna et al. 2019). A local pressure maximum forms near the expanding cavity’s edge, which efficiently traps drifting dust (Pinilla et al. 2018; Gárate et al. 2021, 2023). This evolving pressure bump can enhance the local metallicity and trigger the SI, potentially facilitating late-stage planetesimal formation (Carrera et al. 2017; Ercolano et al. 2017). Similarly, Throop & Bally (2005) pointed out that the external photoevaporation from the nearby massive stars removes gas efficiently, while the dust population remains largely unaffected. This can also increase the disk dust-to-gas ratio.

This paper investigates whether stellar X-ray photoevaporation can produce the conditions required to trigger the SI and form planetesimals. We employed the DustPy code to couple models of viscous gas evolution with dust coagulation, fragmentation, and radial drift, explicitly including photoevaporation to simulate the formation and inside-out movement of a disk cavity with its associated pressure bumps. By applying established SI criteria to these evolving disks, we determined the circumstances under which the photoevaporation-driven dispersal promotes planetesimal formation.

Our paper is organized as follows. The model is introduced in Sect. 2. Our main results and parameter exploration are presented in Sect. 3. The implications and uncertainties are discussed in Sect. 4. The conclusions are drawn in Sect. 5.

2 Method

In this section, we introduce the basics of our model. The prescriptions of disk evolution, stellar X-ray photoevaporation, and planetesimal formation are described in Sects. 2.1, 2.2, and 2.3, respectively, whereas the simulation setup is given in Sect. 2.4.

2.1 Disk model

2.1.1 Gas component

The initial gas disk surface density is adopted from the self-similar solution of Lynden-Bell & Pringle (1974): $${\Sigma }_{\mathrm{g}0}=\frac{M_{\text{disk}}}{2\pi R_{\mathrm{c}}^2}{\left(\frac{r}{R_{\mathrm{c}}}\right)}^{-1}exp(-\frac{r}{R_{\mathrm{c}}}),$$(1) where M_disk is the onset gas disk mass, R_c is the disk characteristic radius, and r is the radial distance from the central star.

Considering a passively irradiated disk with a constant irradiation angle of θ=0.05, the disk midplane temperature can be written as $$T={\left(\frac{\theta L_{\star }}{8\pi r^2{\sigma }_{\mathrm{SB}}}\right)}^{1/4}=220\mathrm{K}{\left(\frac{r}{\mathrm{au}}\right)}^{-1/2},$$(2) where $L_{\star }=4\pi {\sigma }_{\mathrm{SB}}{R}_{\star }^2{T}_{\text{eff}}^4$ is the stellar luminosity, and σ_SB is the Stefan-Boltzmann constant. We chose a pre-main sequence solar mass central star with a stellar radius of R_*=2 R_⊙ and an effective temperature of T_eff=5772 K.

The aspect ratio of the gas disk is therefore expressed as $$h_{\mathrm{g}}=\frac{c_{\mathrm{s}}}{v_{\mathrm{K}}}=\frac{H_{\mathrm{g}}}{r}=0.03{\left(\frac{r}{\mathrm{au}}\right)}^{1/4},$$(3) where c_s is the isothermal gas sound speed, $v_{\mathrm{K}}=\sqrt{G{M}_{\star }/r}$ is the Keplerian velocity, G is the gravitational constant, M_* is the stellar mass, and H_g is the gas disk scale height.

The gas pressure gradient parameter is defined as $$\eta =-\frac{h_{\mathrm{g}}^2}{2}\frac{\partial ln{P}_{\mathrm{g}}}{\partial lnr},$$(4) where P_g is the midplane gas pressure. Another useful dimensionless parameter that quantifies the gas pressure effect can be defined as $$\Pi =\frac{\eta v_{\mathrm{K}}}{c_{\mathrm{s}}},$$(5) where ηv_K is the gas headwind velocity, indicating how much the azimuthal velocity of the gas deviates from the local Keplerian velocity.

The gas surface density is evolved through the viscous advection-diffusion equation (Lynden-Bell & Pringle 1974): $$\frac{\partial }{\partial t}{\Sigma }_{\mathrm{g}}+\frac{1}{r}\frac{\partial }{\partial r}\left(r{\Sigma }_{\mathrm{g}}{v}_{\mathrm{g}}\right)={\dot{\Sigma }}_{\mathrm{PE},\mathrm{g}},$$(6) where v_g is the gas radial velocity and the sink term ${\dot{\Sigma }}_{\text{PE,g}}$ is the gas surface density decay rate driven by the stellar X-ray radiation, the prescription of which is presented in Sect. 2.2.

The radial viscous velocity of the gas reads $$v_{\mathrm{g}}=-\frac{3}{{\Sigma }_{\mathrm{g}}\sqrt{r}}\frac{\partial }{\partial r}\left({\Sigma }_{\mathrm{g}}v\sqrt{r}\right),$$(7) where v=αc_sH_g is the gas turbulent viscosity, and α is the dimensionless viscous parameter (Shakura & Sunyaev 1973).

2.1.2 Dust transport

The advection-diffusion equation for the dust is given by Clarke & Pringle (1988): $$\frac{\partial }{\partial t}{\Sigma }_{\mathrm{d}}+\frac{1}{r}\frac{\partial }{\partial r}[r{\Sigma }_{\mathrm{d}}{v}_{\mathrm{d}}-r{D}_{\mathrm{d}}{\Sigma }_{\mathrm{g}}\frac{\partial }{\partial r}\left(\frac{{\Sigma }_{\mathrm{d}}}{{\Sigma }_{\mathrm{g}}}\right)]={\dot{\Sigma }}_{\mathrm{PE},\mathrm{d}}+{\dot{\Sigma }}_{\mathrm{plt}},$$(8) where Σ_d and v_d are the dust surface density and radial velocity, and ${\dot{\Sigma }}_{\text{PE,}\mathrm{d}}$ and ${\dot{\Sigma }}_{\text{plt}}$ are the dust mass-loss rates due to stellar X-ray photoevaporation and planetesimal formation, respectively. Since the dust coagulation is considered, in practice, the above transport equation is solved for dust with different mass bins. The disk metallicity is defined as Z=Σ_d/Σ_g.

The dust diffusivity taken from Youdin & Lithwick (2007) reads $$D_{\mathrm{d}}=\frac{D_{\mathrm{g}}}{1+{\mathrm{St}}^2}=\frac{{\delta }_{\mathrm{r}}{c}_{\mathrm{s}}{H}_{\mathrm{g}}}{1+{\mathrm{St}}^2},$$(9) where D_g=δ_r c_s H_g is the gas diffusivity and the dimensionless turbulent mixing parameter δ_r represents the turbulent diffusion in the radial direction. We set D_d=0 for most of our simulations. Including the dust diffusivity has a negligible effect on our results, since the timescale for dust to diffuse across its trapping region is much longer than that of planetesimal formation (see discussions in Sect. 4.2).

The aerodynamic size of the dust particles can be quantified by their Stokes number (Whipple 1972, 1973; Weidenschilling 1977): $$\mathrm{S}\mathrm{t}=\{\begin{array}{ll}\frac{\pi }{2}\frac{a_{\mathrm{d}}{\rho }_{\bullet }}{{\mathrm{\Sigma }}_{\mathrm{g}}}& \text{at}{a}_{\mathrm{d}}\le \frac{9}{4}{\lambda }_{\mathrm{m}\mathrm{f}\mathrm{p}},\\ \frac{2\pi }{9}\frac{a_{\mathrm{d}}^2{\rho }_{\bullet }}{{\lambda }_{\mathrm{m}\mathrm{f}\mathrm{p}}{\mathrm{\Sigma }}_{\mathrm{g}}}& \text{at}{a}_{\mathrm{d}}>\frac{9}{4}{\lambda }_{\mathrm{m}\mathrm{f}\mathrm{p}}.\end{array}$$(10) Here a_d is the particle’s physical size, λ_mfp is the mean free path of disk gas, and ρ_• is the particle’s internal density, chosen to be 1.67 g cm⁻³ as the default value.

The dust radial velocity is given by $$v_{\mathrm{d}}=\frac{v_{\mathrm{g}}+2v_{\mathrm{drift},max}\mathrm{St}}{1+{\mathrm{St}}^2},$$(11) with the maximum drift velocity v_drift,max=−ηv_K.

2.1.3 Dust coagulation

In our simulations, the dust coagulation is calculated by solving the Smoluchowski equation (Smoluchowski 1916; Stammler & Birnstiel 2022). Basically, the dust initially grows through pairwise collisions until reaching the radial drift barrier or the fragmentation barrier.

Although our simulations employed the DustPy code (Stammler & Birnstiel 2022), which directly solves the full Smoluchowski equation, we briefly introduce the concepts of the fragmentation and drift barriers to provide physical intuition for the dust growth process. These barriers describe the limiting regimes in which dust growth is halted by either fragmentation or radial drift, and they are commonly used to interpret the size limits of dust in protoplanetary disks.

The dust growth, fragmentation, and radial drift have been investigated in-depth (Safronov 1972; Weidenschilling 1977, 1980; Nakagawa et al. 1986; Dominik & Tielens 1997; Brauer et al. 2008; Blum & Wurm 2008; Birnstiel et al. 2012).

Assuming that the relative velocity between dust is driven by disk turbulence $v_{\text{rel}}\simeq \sqrt{3{\delta }_{\mathrm{t}}\mathrm{St}}{c}_{\mathrm{s}}$ (Ormel & Cuzzi 2007), the growth timescale taken from Birnstiel et al. (2012) reads $${\tau }_{\text{growth}}=\frac{a_{\mathrm{d}}}{{\dot{a}}_{\mathrm{d}}}=\frac{{\rho }_{•}{a}_{\mathrm{d}}}{{\rho }_{\mathrm{d}}{v}_{\text{rel}}}\simeq \frac{1}{Z{\Omega }_{\mathrm{K}}},$$(12) where ρ_d is the dust volume density, and Ω_K=v_K/r is the Keplerian frequency. The dust radial drift timescale adopted from Birnstiel et al. (2012) reads $${\tau }_{\mathrm{drift}}=\frac{r}{v_{\mathrm{d}}}\simeq \frac{1}{\eta \mathrm{St}{\Omega }_{\mathrm{K}}}.$$(13) When the drift timescale is shorter than the growth timescale, the in situ growth cannot proceed. By equating these two timescales, one can derive the particles’ Stokes number in the radial drift barrier as $${\mathrm{St}}_{\mathrm{drift}}=\frac{Z}{2\eta }.$$(14) On the other hand, when the dust particles’ relative velocities exceed a specific threshold, they fragment during collisions rather than perfectly merge. This sets another limit on the size of the particles, referred to as the fragmentation barrier. The fragmentation barrier in the turbulence-driven regime is given by $${\mathrm{St}}_{\mathrm{frag}}=\frac{v_{\mathrm{frag}}^2}{3{\delta }_{\mathrm{t}}{c}_{\mathrm{s}}^2},$$(15) where v_frag is the fragmentation threshold velocity and δ_t is the turbulent mixing parameter. In reality, the Stokes number of the largest particles can be approximated as St=min(St_frag, St_drift).

It is worth noting that the turbulent mixing parameters are defined in three forms: δ_r, δ_t, and δ_z. For MRI-driven disk turbulence, which is approximately isotropic, angular momentum transport is roughly equivalent to mass diffusion (Johansen & Klahr 2005; Zhu et al. 2012; Zhou et al. 2022) (v ≃D_g), so that δ_r ≃ δ_t ≃ δ_z ≃ α. In our default models, we set δ_r=0, δ_t ≃ δ_z ≃ α, and explore the influences of these parameters in Sect. 4.2.

2.2 Stellar X-ray photoevaporation

Intense X-ray radiation from the young, active central star heats the upper layers of the disk, inducing the hot gas and dust to escape as a photoevaporative wind (Alexander et al. 2014; Ercolano & Pascucci 2017). In this work we focus on the stellar X-ray radiation, the influence of which is incorporated as the sink terms in Eqs. (6) and (8).

We adopt the normalized ${\dot{\Sigma }}_{\text{PE,g}}$ profiles from Owen et al. (2012) for primordial disks and disks with inner holes. We also note that alternative prescriptions, such as that of Picogna et al. (2019), yield different radial mass-loss distributions due to differences in their hydrodynamic setups and thermochemical treatments. More recently, Sellek et al. (2024) find that Owen et al. (2012)’s prescription remains in good agreement with their radiation-hydrodynamic simulations. This provides the motivation for our choice of Owen et al. (2012)’s prescription in this work. Although these prescriptions differ in their detailed radial mass-loss profiles, a systematic exploration of how these differences may impact the disk dispersal and dust dynamics is beyond the scope of this study.

Owen et al. (2012) calculated the ${\dot{\Sigma }}_{\text{PE,g}}$ profiles for both primordial disks and disks with inner holes, based on their radiation-hydrodynamical model. $${\dot{M}}_{\mathrm{P}\mathrm{E}}=\{\begin{array}{rll}1.65& \times {10}^{-8}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-0.068}& \\ & \times {\left(\frac{L_{\mathrm{X}}}{L_{\mathrm{X},\odot }}\right)}^{1.14}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}& [\text{\hspace{0.17em}primordial disk\hspace{0.17em}}]\\ 1.27& \times {10}^{-8}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-0.148}& \\ & \times {\left(\frac{L_{\mathrm{X}}}{L_{\mathrm{X},\odot }}\right)}^{1.14}{M}_{\odot }{\mathrm{y}\mathrm{r}}^{-1}& [\text{disk with inner hole,}\\ & & \text{at}r>r_{\text{hole\hspace{0.17em}}}]\end{array}$$(16)

The typical X-ray luminosity of a solar-mass T Tauri star can be expressed as L_{X, ⊙}=2.34 × 10³⁰ erg s⁻¹ (Preibisch et al. 2005). Taking into account the correlation between the X-ray luminosity of T Tauri stars and stellar mass (see Sect. 4.2 in Preibisch et al. 2005), the integrated gas mass-loss rate driven by the stellar X-ray photoevaporation can be written as $${\dot{M}}_{\mathrm{P}\mathrm{E}}=\{\begin{array}{lll}1.65\times {10}^{-8}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{1.57}& M_{\odot }{\mathrm{y}\mathrm{r}}^{-1}& [\text{\hspace{0.17em}primordial disk}],\\ 1.27\times {10}^{-8}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{1.49}& M_{\odot }{\mathrm{y}\mathrm{r}}^{-1}& \text{\hspace{0.17em}[disk with inner hole,}\\ & & \text{at}r>r_{\text{hole}}].\end{array}$$(17)

For the primordial disks, the photoevaporative mass-loss rate ${\dot{\Sigma }}_{\text{PE,g}}$ peaks at ∼2−3 au. After the cavity forms, we set the gas depletion rate at r<r_hole to a constant to ensure rapid gas removal there. Since the inner disk does not participate in planetesimal formation, this treatment does not alter the final mass of planetesimals. We adopted a relatively high X-ray luminosity L_X=6.15 × 10³⁰ erg s⁻¹=2.6 L_{X, ⊙} in our default model and assume that L_X remains constant during the short disk dispersal phase. Our simulations are designed to follow the final dispersal of the disk. Choosing a higher L_X, corresponding to a more X-ray-active young star (see Fig. 3 in Preibisch et al. 2005), allows the photoevaporation-driven clearing phase to complete within ∼1 Myr.

The stellar photoevaporation also influences the dust dynamics. Following Facchini et al. (2016) and Sellek et al. (2020), we assume that particles smaller than the entrained size, a_ent, can be blown away by the wind, where $$a_{\text{ent}}=\frac{v_{\text{th}}{\dot{M}}_{\mathrm{PE}}}{4\pi \mathcal{F}{\rho }_{•}G{M}_{\star }},$$(18) where $v_{\text{th}}=\sqrt{8/\pi }{c}_{\mathrm{s}}$ is the thermal velocity of the gas and $\mathcal{F}=H_{\mathrm{g}}/\sqrt{r^2+H_{\mathrm{g}}^2}$ corresponds to the disk geometric factor. For particles smaller than the entrained size a_ent, we set the dust mass-loss rate: $${\dot{\Sigma }}_{\mathrm{PE},\mathrm{d}}^i=f_{\mathrm{ent}}·Z^i·{\dot{\Sigma }}_{\mathrm{PE},\mathrm{g}},$$(19) where ${\dot{\Sigma }}_{\text{PE,d}}^i$ is the dust mass-loss rate driven by the stellar X-ray photoevaporation within a specific dust mass bin i, f_ent is the entrained mass fraction, and Zⁱ is the metallicity for that mass bin. Following Sellek et al. (2020), we set f_ent=1 in our simulations.

2.3 Streaming instability criterion and planetesimal formation

Several literature studies have investigated the critical condition for triggering the SI. The classical criterion for the onset of the SI is proposed by Youdin & Goodman (2005) through a linear stability analysis. They found that the SI criterion can be fulfilled when the midplane volume density ratio of dust-to-gas approaches unity. Their result laid the foundation for subsequent theoretical and numerical studies. Carrera et al. (2015) explored this issue by performing 2D numerical simulations and obtained the metallicity threshold Z_crit as a function of particles’ Stokes number. Later, Yang et al. (2017) reduced this clumping threshold with an improved numerical method. Li & Youdin (2021) further quantified how Z_crit(St) varied with simulation setups. They found that subsolar metallicity can already induce the dust clumping at St∼0.1-01, lower than previous results. However, their numerical study was only 2D and did not take into account disk turbulence. More recently, Lim et al. (2024) advanced the understanding of the SI condition by performing 3D stratified shearing box simulations that included both particles’ self-gravity and externally forced turbulence. Their work provides a more accurate threshold for evaluating the planetesimal formation in global disk models.

For the purpose of generality, here we list three SI criteria, ranging from the early linear analysis of Youdin & Goodman (2005) to the most updated numerical simulations of Lim et al. (2024). The influence of adopting different SI criteria on our study is discussed in Sect. 3.1.

To the zero order, the analytical SI criterion can be expressed as the midplane volume density ratio of dust-to-gas, where $${ϵ}_{\text{crit}}\equiv {\rho }_{\mathrm{d}}/{\rho }_{\mathrm{g}}=1.$$(20)

By explicitly quantifying how turbulence affects the SI-induced clumping, the new derived criteria from Lim et al. (2024) bridge the gap between idealized SI theory and turbulent disk environments. The ε_crit criterion provided by Lim et al. (2024) reads $${log}_{10}{ϵ}_{\text{crit}}=0.42{\left({log}_{10}\mathrm{St}\right)}^2+0.72\left({log}_{10}\mathrm{St}\right)+0.37.$$(21) Meanwhile, Lim et al. (2024) also provides the Z_crit criterion: $$\begin{array}{rl}{\log }_{10}{Z}_{\text{crit}}=& 0.15{\left({\log }_{10}{\delta }_{\mathrm{t}}\right)}^2-0.24({\log }_{10}\mathrm{S}\mathrm{t}{\log }_{10}{\delta }_{\mathrm{t}})\\ & -1.48({\log }_{10}\mathrm{S}\mathrm{t})+1.18\left({\log }_{10}{\delta }_{\mathrm{t}}\right)\end{array}$$(22) The above two criteria hold for particles with 0.01 ≤ St ≤ 0.1 and a disk turbulent level of 10⁻⁴ ≤ δ_t ≤ 10⁻³. We also use YG05, L24E, and L24Z to represent Eqs. (20), (21), and (22) hereafter.

It is worth noting that L24Z (Eq. (22)) is derived at Π=0.05. Physically, a more decisive factor that sets the SI clumping should be Z/Π in the above equation (Sekiya & Onishi 2018). In our simulations that include the stellar X-ray photoevaporation, however, the Π value varies significantly across the disk and crosses zero at the outer cavity edge. This variability makes it impractical to apply a correction to the criterion directly. We therefore adopt the ε_crit criterion of L24E (Eq. (21)), which shares the same physical basis as the Z_crit criterion of L24Z (Eq. (22)), and subsequently calculate the midplane density ratio in our 1D code to incorporate the effects of Π. A key advantage of L24E (Eq. (21)) is its ability to capture turbulence variations from the SI (specifically through Π) via the self-consistent modification of the dust scale height. The specific method is detailed below.

In the absence of the SI-induced dust concentration, the dust scale height adopted from Dubrulle et al. (1995) can be written as $$H_{\mathrm{d}}=H_{\mathrm{g}}\sqrt{\frac{{\delta }_{\mathrm{z}}}{{\delta }_{\mathrm{z}}+\mathrm{St}}},$$(23) where δ_z is the vertical turbulent mixing parameter. SI can act as an additional turbulent source to affect the vertical dust distribution. Taking this into account, Carrera et al. (2025) and Eriksson et al. (2025) derive the effective dust scale height as $$H_{\mathrm{d}}=H_{\mathrm{g}}\sqrt{\frac{1+\mathrm{St}}{1+\mathrm{St}+ϵ}}\sqrt{{\left(\frac{\Pi }{5}\right)}^2+\frac{{\delta }_{\mathrm{z}}}{{\delta }_{\mathrm{z}}+\mathrm{St}}}.$$(24) Our simulation is 1D and only calculates the vertically integrated gas and dust quantities. To derive the volume density, the vertical distribution of the dust needs to be intrinsically assumed. Despite that the dust vertical distribution might not be Gaussian-like owing to the SI effect (Lim et al. 2024), we still assume that the dust midplane density approximately follows ρ_d ∝ Σ_d/H_d, and therefore, $$ϵ\equiv \frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}=Z\frac{H_{\mathrm{g}}}{H_{\mathrm{d}}}.$$(25)

Let $\xi \equiv \frac{Z}{\sqrt{1+\mathrm{St}}\sqrt{{\left(\frac{\Pi }{5}\right)}^2+\frac{{\delta }_{\mathrm{z}}}{{\delta }_{\mathrm{z}}+\mathrm{St}}}}$, we can obtain $$ϵ=\frac{{\xi }^2+\sqrt{{\xi }^4+4{\xi }^2(1+\mathrm{St})}}{2}.$$(26)

We note that the aforementioned SI criteria are all derived from simulations based on single-size particles (see Bai & Stone 2010; Zhu & Yang 2021 for multispecies). In this study we consider the drift and coagulation of dust particles with a wide size range. We thus use the mass-weighted mean Stokes number to represent the bulk dust population, the formula of which can be written as $${\mathrm{St}}_{\text{mean}}=\frac{{\sum }_i{\mathrm{St}}^i·{\Sigma }_{\mathrm{d}}^i}{{\sum }_i{\Sigma }_{\mathrm{d}}^i},$$(27) where Stⁱ and ${\Sigma }_{\mathrm{d}}^i$ are the corresponding particles’ Stokes number and mass density in the mass bin i. The Stokes number of the dust population can therefore be specified at any given t and r. We adopt St_mean for calculating ε_crit and H_d.

Once ε ≥ ε_crit is satisfied, we transform the dust into planetesimals via $${\dot{\Sigma }}_{\mathrm{plt}}=-\zeta \frac{{\Sigma }_{\mathrm{d}}}{t_{\mathrm{set}}}=-\zeta {\Sigma }_{\mathrm{d}}\frac{{\Omega }_{\mathrm{K}}}{1+{\mathrm{St}}^{-1}},$$(28) where t_set=[Ω_K /(1+St⁻¹)]⁻¹ is the settling time, which corresponds to the timescale on the formation of the SI filaments (Yang et al. 2017). We adopt this formula of t_set to eliminate interference of particles with extremely large Stokes numbers, such as the particles within the inner cavity (see Fig. 1b, Sect. 3.1).

The planetesimal conversion efficiency, ζ, represents the fraction of dust converted into planetesimals per settling time. Numerically, if ζ is too large, the entire dust mass within a radial bin would be depleted within a single timestep, leading to numerical instability. Physically, once the SI condition is met, dust trapped in a pressure bump is expected to rapidly transform into planetesimals within dozens of settling times (t_set). Given that t_set is significantly shorter than the timescale for the movement of the photoevaporative cavity (∼10³−10⁵ yr), planetesimal formation should be efficient regardless of the adoption of ζ values across a broad parameter range. Based on the simulations of Simon et al. (2016) and following Schoonenberg et al. (2018), we set ζ=0.1. For completeness, we test several planetesimal formation prescriptions. Our fiducial model adopts the rate modified from Schoonenberg et al. (2018). We further examine the formulations by Drążkowska et al. (2016), which consider particles with St>0.01 and an efficiency parameter ζ=0.01, as well as the smooth transition function introduced by Miller et al. (2021) to determine the local planetesimal formation fraction. The resulting differences in the total planetesimal mass are minor, suggesting that our results are robust against the specific choice of the planetesimal formation formula.

We also note that although the SI criterion of L24E (Eq. (21)) is derived for particles with Stokes numbers in the range of ∼0.01−0.1, in our study the dominated dust population involved in planetesimal formation naturally falls within a similar Stokes number range constrained by the drift and fragmentation barriers. The particles’ mean Stokes number (St_mean) reaches only approximately a few tenths even during the late stages of disk evolution (see Fig. 1b, Sect. 3.1). Practically, we extrapolate the thresholds beyond their original Stokes number range. We also test this approach and find that the difference is minor when restricting the Stokes number to ∼0.01−0.1 compared to using the full-size domain.

Fig. 1 Time and radial distance evolution of different quantities in the fiducial model. Panel a: gas (dashed lines) and dust (solid lines) surface density. Panel b: Stokes number of the maximum-size dust particles (solid lines) and the mean Stokes number of the dust particles (dashed lines) at each radial bin. Maximum-size refers to the size of the particle whose surface density is at the peak of the mass bins hereafter. Panel c: midplane density ratio ε (solid lines) and metallicity Z (dashed lines).

2.4 Simulation setup

We used the DustPy (v1.0.5) code (Stammler & Birnstiel 2022) to study the radial evolution of the gas and the dust in protoplanetary disks, with the effects of viscous gas evolution, stellar X-ray photoevaporation, dust coagulation, fragmentation, and radial drift. The radial grid extended from 1–300 au with 256 logarithmically spaced grid cells. We additionally implemented the local refinement of the grids in the planetesimal formation region driven by stellar X-ray photoevaporation, spanning 5–20 au in the typical simulation setups. We verified that the final mass of planetesimals converged with this treatment and a higher radial grid resolution. We adopted 120 bins for dust mass, logarithmically from 10⁻¹² to 10⁵ g. The inner boundary conditions for the dust and gas densities were set to be constant gradients, whereas the outer boundary conditions were set to be floor values.

There are six key parameters in our model: the initial metallicity Z₀, the stellar X-ray luminosity L_X, the gas viscous parameter α, the initial gas disk mass M_disk, the disk characteristic disk radius R_c, and the dust fragmentation threshold velocity v_frag. The default parameter values are listed in Table 1.

3 Results

In this section, we describe the results of the fiducial run in Sect. 3.1 and explore different model parameters from Sects. 3.2 to 3.6. Resultant planetesimal masses and dust-to-planetesimal conversion efficiencies are listed in Table 2, and the radial mass distributions of planetesimals are presented in Fig. 5.

3.1 Fiducial model

For the fiducial run (run_fid), we chose initial disk mass M_disk=0.05 M_⊙, characteristic disk radius R_c=60 au, viscous parameter α=10⁻⁴, initial metallicity Z₀=0.01, fragmentation threshold velocity v_frag=5 m s⁻¹, and stellar X-ray luminosity L_X=2.6 L_{X, ⊙}. The SI criterion is adopted from L24E (Eq. (21)).

Fig. 1 a shows the dust and gas disk evolution at different snapshots in run_fid. At t=0.35 Myr, an inner cavity opens at around 3 au. The outer edge of this inner cavity expands outward rapidly, reaching ∼10 au at t=0.45 Myr, and ∼100 au at t=0.90 Myr, respectively. Due to the dust radial drift and gas disk inside-out clearing, the dust gets accumulated and forms a ring-like structure at the edge of the inner disk cavity.

This can also be clearly seen in Fig. 2, which presents the synthetic dust continuum images at the 880 μm wavelength at three different epochs. At the early phase of t=0.25 Myr, the disk is driven by viscous evolution, and the dust is distributed relatively smoothly throughout the disk (see Fig. 2a). Fig. 2b shows the onset of photoevaporation-driven clearing, whereas in Fig. 2c, the inner cavity expands further out, with a strong dust emission at its outer edge. At the cavity’s edge, rapid gas clearing yields a local pressure bump, providing a natural site for dust concentration. The images were generated with RADMC-3D based on the one-dimensional DustPy outputs mapped onto two-dimensional axisymmetric grids. We adopted 38 logarithmically spaced dust-size bins between 5 μm and 17 cm, and assumed that the disk is passively heated by stellar irradiation. Dust temperatures were computed by Monte Carlo radiative transfer, and synthetic images were convolved to 880 μm, using the DSHARP opacity model (Birnstiel et al. 2018).

Fig. 1b demonstrates the Stokes number of the maximum-size dust particles (solid lines) and the mean Stokes number of the dust population (dashed lines). Overall, the St (a_max) and St_mean values are higher near the outer edge of cavity than in the residual disk regions. This is because dust drifts slower near the pressure bump, enhancing the local coagulation. We also obtain that St_mean near the cavity’s edge is ∼0.05 at t=0.45 Myr, ∼0.1 at t=0.60 Myr, and ∼0.2 at t=0.90 Myr. Thus, St_mean remains largely below 0.1 within 30 au, which is the most relevant for planetesimal formation. This size is within the range where L24E (Eq. (21)) is valid. We also note that the sharp increase in St_mean just inside the cavity is caused by the floor values of the surface densities set in DustPy. These anomalous Stokes number values, as discussed in Sect. 2.3, can lead to the artificial formation of planetesimals in dust-depleted regions. Our adopted formula for t_set (Eq. (28)) is designed to prevent this unphysical outcome.

Fig. 1c illustrates the ratios of dust-to-gas volume density ε and surface density Z. Due to dust trapping at the pressure bump, the disk metallicity Z gets largely elevated and exceeds 2−4%, whereas ε reaches the order of unity at the edge of the inner cavity as it sweeps outward.

We adopt the SI criterion from L24E (Eq. (21)). In the end, planetesimals with a total mass of 31.4 M_⊕ are generated through this photoevaporation gas clearing process. The dust-to-planetesimal conversion efficiency is 20.4%, which indicates that 20.4% of the initial dust mass has been converted into planetesimals. The influences of alternative criteria of YG05 (Eq. (20)) and L24Z (Eq. (22)) on planetesimal formation are also tested. The total masses of resultant planetesimals show minor differences among all these criteria (≲ 10%). Thus, the ε_crit criterion of L24E is used as the default in all following simulations.

Fig. 3a shows the evolution of the total dust mass (black line) and the planetesimal mass (red line) in run_fid. The blue bars represent the mass distribution of dust across different size bins within the pressure bump. The rapid rise of the planetesimal mass (red line) indicates that trapped dust particles are efficiently converted into planetesimals. This high formation efficiency results from the midplane volume density ratio ε being sufficiently high to trigger the SI at the cavity edge. Furthermore, after cavity opening, the dust population within the pressure bump is predominantly composed of particles with Stokes numbers of 0.01−0.1. The contribution of the dust particles within the pressure bump to the final planetesimal mass depends solely on the mass fraction of the trapped dust particles.

To conclude, run_fid shows that stellar X-ray photoevaporation opens a cavity at orbital radii of a few au. As the inner cavity sweeps inside-out, the dust particles accumulating at the cavity’s edge can trigger the SI and form planetesimals efficiently. In total, 30.5 M_⊕ of planetesimals are generated at r=5−30 au within 1 Myr after the onset of stellar X-ray photoevaporation.

Fig. 2 Synthetic 880 μm continuum images from the fiducial model. Panel a: t=0.25 Myr, pre-cavity stage; Panel b: t=0.35 Myr, onset of photoevaporative cavity opening; Panel c: t=0.45 Myr, post-cavity stage.

3.2 Metallicity

We investigated the influence of disk metallicity on planetesimal formation by increasing Z₀ from 0.01 in run_fid to 0.02 in run_metal. The other parameters are kept identical as run_fid.

Fig. 3b shows total dust mass and planetesimal mass in run_metal. The initial dust mass is ∼300 M_⊕, whereas 50 M_⊕ remains when the photoevaporation opens a cavity at t=0.37 Myr. For comparison, in run_fid, the initial dust mass is ∼150 M_⊕, whereas 37 M_⊕ remains at t=0.37 Myr (see Fig. 3a).

We find that run_fid yields a final planetesimal mass of 31.4 M_⊕, whereas in run_metal the final planetesimal mass reaches 48.2 M_⊕. Clearly, doubling the initial metallicity indeed results in an increase in the final planetesimal mass. But the final planetesimal mass in run_metal is not twice as high as in run_fid. This is because the dust in the run_metal experiences more rapid growth (see Eq. (12)) compared to that in run_fid. In the early phase, the grown particles also get accreted by the central star at a faster rate. As a consequence, in both cases, the masses of retained dust at the onset of photoevaporative cavity are overall similar.

3.3 Stellar X-ray luminosity

Meanwhile, we also tested the effect of the stellar X-ray luminosity such that L_X was set to be 1.0 L_{X, ⊙} in run_lumi rather than 2.6 L_{X, ⊙} in run_fid. The other parameters are the same as run_fid.

Table 2 shows that run_fid yields a final planetesimal mass of 31.4 M_⊕ with an overall dust-to-planetesimal conversion efficiency of 20.4%, whereas run_lumi achieves a final planetesimal mass of only 1 M_⊕ with a low conversion efficiency of 0.6%. As can be seen from Fig. 3c, only 3.8 M_⊕ of dust is left when the photoevaporation opens a cavity at t=1.2 Myr in run_lumi.

The low stellar X-ray luminosity model results in a significant drop in the planetesimal mass and conversion efficiency. The strength of L_X is crucial to determine the time for the photoevaporation-driven cavity formation. The disk around a star with a low stellar X-ray luminosity results in the cavityopening at a later time. By this time, a substantial amount of dust has already been accreted to the star. As a consequence, only a limited amount of dust remains available for the planetesimal formation. Therefore, we suggest that a relatively high L_X is required to induce the generation of planetesimals at the gas clearing phase.

Fig. 3 Evolution of dust mass (black lines), planetesimal mass (red lines), and dust mass fraction within the pressure bump in three Stokes number bins (light blue for St<10−2, blue for 10−2 ≤ St ≤ 10−1, and deep blue for St>10−1), across four evolutionary stages: (i) pre-cavity, (ii) cavity at 10 au, (iii) cavity at 20 au, and (iv) cavity at 100 au. These panels show how initial conditions affect dust retention and planetesimal formation: Panel a: fiducial model (Z=0.01, LX=2.6 LX, ⊙, Mdisk=0.05 M⊙, Rc=60 au, α=10−4) serves as the baseline for comparison, showing balanced dust evolution and planetesimal production across the disk. Panel b: higher metallicity (Z=0.02) facilitates the retention of dust mass, thereby promoting planetesimal formation. Panel c: decreased X-ray luminosity (from 2.6 LX, ⊙ to 1.0 LX, ⊙) diminishes dust mass and consequently impedes planetesimal formation. Panel d: higher viscosity parameters (α=10−3, δt=10−4) induce gas replenishment from the outer disk, delay cavity opening, reduce the overall dust reservoir, and thus hinder planetesimal formation. Panel e: a more massive and extended disk (Mdisk=0.1 M⊙, Rc=120 au) prolongs dust retention, thereby extending the duration of planetesimal formation. Panel f: a less massive but extended disk (Mdisk=0.05 M⊙, Rc=120 au) lowers the inner-disk dust density, causing photoevaporation to open a cavity earlier. This makes more dust retained at the onset of photoevaporative cavity, which also provides more building blocks for planetesimal formation.

3.4 Viscosity

Disk viscosity, parameterized by α, drives the global evolution of the protoplanetary disk. We focus on isolating the role of viscosity α in global disk evolution and that of the turbulent mixing parameter δ_t in dust coagulation. In other words, we treat α and δ_t differently. In run_alpha we only varied α from 10⁻⁴ to 10⁻³ whereas the turbulent mixing parameter δ_t was still kept the same. Thus, the particles’ Stokes number set by the fragmentation limit is the same as in run_fid.

Fig. 3d reveals that at a high α of 10⁻³, owing to more efficient angular momentum transport, the gas from the outer disk region replenishes the inner disk region at a fast pace. The cavity is opened at a later time of 0.45 Myr in run_alpha. As a result, a smaller amount of dust is left by this time and eventually leads to less planetesimal formation. Table 2 shows that the planetesimal mass decreases from 31.4 M_⊕ in run_fid to 14.6 M_⊕ in run_alpha.

Fig. 4 Distribution of dust as a function of radial distance and particle size at t=0.30 Myr for run_fid (left) and run_frag (right). The color bar represents the surface density of dust per logarithmic mass interval. The green line is the drift limit (Eq. (14)), whereas the blue line is the fragmentation limit (Eq. (15)). The four white lines correspond to St=[10−3, 10−2, 10−1, 100].

3.5 Disk radius

The disk size is also a key factor that affects the gas and dust evolution and planetesimal formation. In run_fid, the characteristic disk radius R_c was set to be 60 au. Here we performed two additional sets of simulations to explore the influence of disk size. In run_dsize1, we kept Σ_g0 at 1 au the same as that of run_fid and increased R_c to 120 au. In such a case the initial disk mass increased to 0.1 M_⊙. On the other hand, in run_dsize2, the total disk mass was still kept the same as 0.05 M_⊙, whereas R_c was increased to 120 au. This configuration leads to an overall lower Σ_g 0 in the inner disk region.

Figs. 3e and f depict the evolution of dust masses and planetesimal masses in these two cases. The disk is more massive, and a higher proportion of dust mass resides in the outer disk region in run_dsize1 compared to run_fid. Since the gas density remains approximately the same in the inner disk region for run_fid and run_dsize1, the photoevaporative cavity forms at a similar time. In run_dsize1, the residual dust mass reaches ∼144 M_⊕ at this cavity-opening time (see Fig. 3e). This leads to a higher planetesimal mass and conversion efficiency in the end.

On the other hand, run_dsize2 refers to a circumstance in which the gas mass is more readily distributed at further out disk regions. For run_dsize2, the gas density in the inner disk region is lower, and the photoevaporation opens a cavity at an earlier time of 0.18 Myr compared to run_fid (see Fig. 3f). This makes more dust retained at the onset of photoevaporative cavity, which also provides more building blocks for planetesimal formation.

As can be seen from Table 2, the conversion efficiencies in the above two models are quite similar, but both are much higher than run_fid. This again indicates that regardless of the detailed treatment of Σ_g or M_disk, a disk with a larger size promotes the dust retention and thus subsequent planetesimal formation.

3.6 Dust fragmentation threshold velocity

To explore how the photoevaporation-driven planetesimal formation depends on the fragmentation-limited dust growth, we performed a simulation run_frag with a lower fragmentation threshold velocity of v_frag=1 m s⁻¹. The other parameters are the same as run_fid.

Fig. 4 shows the distribution of dust as a function of radial distance and particle size at t=0.30 Myr for run_fid (left) and run_frag (right). It is evident that the particles’ Stokes number in run_fid is more than one order of magnitude larger than that in run_frag. Because the growth is limited by fragmentation, the largest particles reach a lower Stokes number when the fragmentation threshold is lower. Such small particles are more strongly coupled to the disk gas and hence drift more slowly. Consequently, they are less efficient at concentrating at the cavity’s edge, preventing the SI clumping and subsequent planetesimal formation. As can be seen from Table 2, no planetesimal can form during the gas clearing phase in run_frag.

Lim et al. (2025b) and Lim et al. (2025a) conducted 2D and 3D simulations of pure SI without external turbulence, respectively. The SI criteria presented in these two studies suggest lower thresholds for the onset of the SI. However, it is worth noting that their simulations indicate that ε needs to exceed 2 to trigger the SI for particles with St ≤ 10⁻³. Hence, even adopting a loose SI criterion from Lim et al. (2025a) (their Eq. (11)) without considering any turbulence, planetesimals are still unlikely to form at a relatively low fragmentation velocity.

Zhao et al. (2025) also concludes that v_frag=2 m s⁻¹ is the minimum value that allows planetesimal formation, which is consistent with our results.

To facilitate a direct comparison among different models, Fig. 5 presents the radial mass distributions of planetesimals. To conclude, our results indicate that planetesimal formation is more efficient in larger disks with higher metallicities, lower viscosities, higher dust fragmentation threshold velocities, and/or around stars with higher X-ray luminosities.

Fig. 5 Comparison of the radial mass distributions of planetesimals at different disk regions. The curves show the results for models with varying parameters: run_fid (black), run_metal (blue), run_lumi (orange), run_alpha (purple), run_dsize1 (red), and run_dsize2 (cyan).

4 Discussions

We compare our work with the literature studies in Sect. 4.1. The influences of turbulent mixing parameters and pressure bumps produced by other mechanisms are presented in Sect. 4.2 and Sect. 4.3. Finally, the caveats in our model are discussed in Sect. 4.4.

4.1 Comparison with literature studies

We compare our results with those of Carrera et al. (2017) and Ercolano et al. (2017), who also explored planetesimal formation via the SI driven by the stellar X-ray photoevaporation. First, our treatment of dust evolution differs from theirs. Whereas Carrera et al. (2017) used a custom 1D viscous disk code with simplified dust prescriptions from Birnstiel et al. (2011), and Ercolano et al. (2017) employed the two-population model of Birnstiel et al. (2012) in their 1D code, we used the open-source DustPy code, which incorporates full multi-size dust evolution, including coagulation, fragmentation, and radial drift. The other difference is the treatment of stellar X-ray photoevaporation. Carrera et al. (2017) adopted the model from Gorti et al. (2015), who considered FUV-driven disk dispersal. This process triggers strong radial pressure gradients and enhances the dust pile-ups in the early phase of the disk lifetime. Ercolano et al. (2017), by contrast, used both X-ray and EUV-driven photoevaporation prescriptions derived from Owen et al. (2010, 2012)’s hydrodynamic simulations. Consequently, in their works, disk gas was removed more gradually, and photoevaporation primarily affected the inner region of the disk. While also adopting the stellar X-ray photoevaporation prescription from Owen et al. (2012), we additionally incorporate dust removal via photoevaporation, following the studies of Facchini et al. (2016) and Sellek et al. (2020).

There are two improvements in our study: the dust scale height calculation and the SI criterion. Whereas the previous studies relied on simplified estimates for the midplane density, we implement the correction from Carrera et al. (2025) and Eriksson et al. (2025) (our Eq. (24)), which self-consistently accounts for local turbulence and particle-size-dependent settling. This allows us to track dust evolution more precisely and significantly improves the accuracy of midplane density predictions. On the other hand, both Carrera et al. (2017) and Ercolano et al. (2017) adopted the classical threshold from the 2D unstratified simulations of Carrera et al. (2015). In this study, we instead adopt the latest criterion from Lim et al. (2024), which is based on 3D shearing box simulations and defines a critical midplane density ratio ε_crit. This criterion considers the external turbulence, reflecting a more realistic condition for dust-gas interaction in vertically stratified disks.

Regarding the planetesimal formation, Carrera et al. (2017) reports up to 60 M_⊕ of planetesimals form beyond the 100 au region and ∼16 M_⊕ within the 100 au region. This high efficiency is a consequence of their adoption of strong FUV photoevaporation, which accelerates disk dispersal in the outer disk region and thus enhances dust concentration there. Conversely, Ercolano et al. (2017) finds planetesimal formation is less likely for the X-ray and EUV photoevaporation. Their model produces only 3 M_⊕ under optimal conditions, due to rapid dust drift before gap formation and absence of strong pressure traps. Our results strike a middle ground. We demonstrate that stellar X-ray photoevaporation can still trigger localized SI when vertical dust settling is explicitly modeled. Because X-ray–driven winds have higher temperatures and thus smaller gravitational radii, the photoevaporative cavity opens at a smaller physical radius, allowing a larger fraction of solids to remain inside ∼100 au and thereby promoting planetesimal formation within ∼100 au. In our fiducial model, the total planetesimal mass can attain 30.5 M_⊕ within the 30 au disk region. By adopting a more realistic SI criterion and a sophisticated dust evolution treatment, our study indicates that stellar X-ray photoevaporation can also be robust for SI-triggered planetesimal formation.

4.2 Turbulent mixing parameters

Several turbulent mixing parameters are used in this work: δ_r, δ_t and δ_z. Since disk turbulence can be nonisotropic, for example, driven by the vertical shear instability (VSI) (Nelson et al. 2013) or nonideal MHD (Yang et al. 2018), these parameters are not necessarily equal. As shown in Table 3, we conducted several simulations to investigate the influences of varying these parameters.

We performed a simulation run_rad, considering the effect of turbulent radial diffusion on dust and other parameters kept identical as run_fid. The dust diffusion timescale in run_rad is given by $$t_{\mathrm{diff}}=\frac{{\left(\Delta r\right)}^2}{D_{\mathrm{d}}}\simeq 1.5\times {10}^3{f}^2{\left(\frac{{\delta }_{\mathrm{r}}}{{10}^{-4}}\right)}^{-1}\left(\frac{P_{\mathrm{r}}}{1\mathrm{yr}}\right)\mathrm{yr},$$(29) where Δ r ∼fH_g is the pressure bump width at the onset of photoevaporative cavity, f is an order-unity factor, and P_r is the Keplerian period at the radial distance r.

The photoevaporative cavity expansion timescale can be estimated as $$t_{\mathrm{cav}}=\frac{\Delta r}{v_{\mathrm{cav}}}=\frac{\Delta r}{r}\frac{{\Sigma }_{\mathrm{cav}}}{{\dot{\Sigma }}_{\mathrm{PE},\mathrm{g}}}\simeq {10}^3{f}^2{\left(\frac{\alpha }{{10}^{-4}}\right)}^{-1}\left(\frac{P_{\mathrm{r}}}{1\mathrm{yr}}\right)\mathrm{yr},$$(30) where $v_{\text{cav}}=r{\dot{\Sigma }}_{\text{PE,g}}/{\Sigma }_{\text{cav}}$ is the cavity expansion velocity, we emphasize that v_cav does not correspond to the physical radial velocity of the gas; rather, it represents the expansion speed (for example, the phase velocity) of the photoevaporative cavity edge. Σ_cav ≃Ṁ_acc/(3πv) denotes the critical gas surface density at which viscous accretion balances photoevaporative mass loss, and ${\dot{M}}_{\text{acc}}\simeq 2\pi r\Delta r{\dot{\Sigma }}_{\text{PE,g}}$ corresponds to the critical gas accretion rate.

We note that, in estimating t_cav, we have assumed ${\dot{M}}_{\text{acc}}\simeq 2\pi r\Delta r{\dot{\Sigma }}_{\text{PE,}\mathrm{g}}$ for simplicity. In practice, however, the gas accretion rate after the onset of cavity opening is typically smaller than the photoevaporative mass-loss rate. This implies that our estimate of t_cav should be regarded as an upper limit. Therefore, we provide only an estimate that t_cav ≲ t_diff. Although the exact reduction factor is difficult to quantify, this effect further reinforces our conclusion that radial diffusion cannot regulate the dust concentration at the cavity edge.

The dust particles settle and are subsequently converted into planetesimals on a timescale of: $$t_{\mathrm{plt}}=\frac{1}{\zeta \mathrm{St}{\Omega }_{\mathrm{K}}}\simeq 16{\left(\frac{\zeta }{0.1}\right)}^{-1}{\left(\frac{\mathrm{St}}{0.1}\right)}^{-1}\left(\frac{P_{\mathrm{r}}}{1\mathrm{yr}}\right)\mathrm{yr},$$(31) where St is the Stokes number of the maximum-size dust particles St (a_max), which follows St(a_max) ≃ 0.1(r/au)^1/2 in the fiducial model.

It can be noted that, for the particles with Stokes numbers of 0.01−0.1, t_plt ≪ t_cav<t_diff in run_rad. Once the SI criterion is met, dust is converted into planetesimals on a timescale much shorter than both the outward expansion of the photoevaporative cavity edge and the dust radial diffusion timescale. In other words, planetesimal formation is completed before the cavity edge recedes and before diffusion can significantly redistribute the dust. Consequently, we find no substantial difference in the total planetesimal mass formed between run_rad and run_fid (see Table 3).

We set δ_r=0, since radial diffusion does not regulate the peak dust-to-gas ratio in the low δ_r and high L_X regime we focus on here. In this regime, photoevaporation drives the cavity edge outward on a timescale shorter than that required for dust radial diffusion to establish a local drift-diffusion equilibrium at the pressure maximum. More importantly, both t_cav and t_diff are much longer than t_plt. Once the SI is triggered, dust is rapidly converted into planetesimals, which naturally limits the peak dust-to-gas ratio.

Following our analysis of the impact of δ_r, we focus on the turbulent mixing parameter δ_t. Furthermore, to investigate the influence of δ_t, we conducted a control simulation, run_turb, based on run_alpha. This control run maintained all other parameters the same as run_alpha, but adjusted δ_t to 10⁻³ and v_frag to 15 m s⁻¹. This ensures that the dust reaches a similar Stokes number at the fragmentation limit in two simulations.

We can therefore isolate the effect of δ_t. Through this comparison, we find that an increase in δ_t leads to a reduction of final planetesimal mass (from 14.6 M_⊕ to 0.7 M_⊕). This phenomenon can be attributed to the fact that the relative velocity between particles increases with δ_t, and dust grows much faster in this circumstance (see Eq. (12)). Consequently, the dust in the inner disk reaches the fragmentation limit at an earlier stage, leading to a less efficient retention of dust mass before the onset of photoevaporative cavity.

Similar to δ_t, δ_z also indirectly affects the planetesimal formation process by altering the dust growth timescale. We also established a control simulation run_vert based on run_fid. We reduced the value of δ_z from 10⁻⁴ in run_fid to 10⁻⁶ in run_vert. This reduction decreases the dust scale height, which both increases the midplane dust-to-gas and shortens the growth timescale. However, ε is not high enough to trigger the SI at early times, and fast growth leads to rapid depletion of dust mass. In the end we find that less dust is retained by the time that cavity opens, leading to a decrease in the final planetesimal mass from 31.4 M_⊕ to 20.9 M_⊕.

4.3 Local pressure bump

Apart from the inner cavity generated by stellar X-ray photoevaporation, the pressure bumps associated with substructures can also be produced by other mechanisms, such as gaps opened by massive planets (Kobayashi et al. 2012; Zhu et al. 2012; Gonzalez et al. 2015; Eriksson et al. 2020), at the inner edge of dead zones (Varnière & Tagger 2006; Lyra et al. 2009; Dzyurkevich et al. 2010), by zonal flows arising from hydrodynamic instabilities (Lehmann & Lin 2025) or magnetorotational turbulences (Johansen et al. 2009a), and in the ice line regions (Ida & Guillot 2016; Drążkowska & Alibert 2017; Schoonenberg & Ormel 2017). The pressure bumps act as traps that halt the inward drift and concentrate dust particles. Under favorable conditions, this process can trigger the SI and lead to planetesimal formation.

To further investigate the role of pressure bumps that may be produced by other mechanisms, we performed an additional simulation in which gaps were artificially imposed on the disk. The gap profile was adopted from Dullemond et al. (2018): $$\alpha \left(r\right)={\alpha }_0/F\left(r\right),$$(32) where α₀ is the disk viscosity, and F(r) is given by $$F\left(r\right)=exp[-Aexp(-{(r-r_0)}^2/\left(2w^2\right))],$$(33) where A is the gap amplitude, r₀ is the gap location, and w is the gap width.

Fig. 6 shows three snapshots of the dust continuum images at 880 μm with two artificial bumps at 30 au and 60 au, where A=0.4 and w=0.1 r₀. Other parameters are the same as run_fid. At t=0.38 Myr, the two wide gaps block the inward drifting of dust from the outer disk region before the onset of photoevaporative cavity (see Fig. 6a). We can see that the photoevaporative cavity forms and progressively moves outward. At t=0.60 Myr the cavity moves to 18 au (see Fig. 6b), whereas at t=0.80 Myr it merges with the inner bump at 30 au (see Fig. 6c).

The bumps can also be sites for dust trapping and planetesimal formation. Consequently, the total planetesimal mass increases from 31.4 M_⊕ to 33.9 M_⊕. Meanwhile, ∼7.6 M_⊕ mass of planetesimals form in the bumps, leading to a slight reduction in the planetesimal mass formed at the outer edge of the photoevaporative cavity, from ∼30.5 M_⊕ to ∼26.3 M_⊕. We also note that as the photoevaporative cavity expands progressively outward, it merges with the existing bumps.

We compute the optical depth τ_v following: $${\tau }_v={\sum }_i({\kappa }_{v,i}^{\mathrm{abs}}+{\kappa }_{v,i}^{\mathrm{sca},\mathrm{eff}}){\Sigma }_{\mathrm{d}}^i,$$(34) where ${\kappa }_{v,i}^{\mathrm{abs}}$ and ${\kappa }_{v,i}^{\mathrm{sca},\text{eff}}=(1-{\mathrm{g}}_{v,\mathrm{i}}){\kappa }_{v,\mathrm{i}}^{\mathrm{sca}}$ are the absorption and effective scattering opacities of dust species i, and ${\Sigma }_{\mathrm{d}}^i$ is the dust surface density of dust species i. The values of ${\kappa }_{v,i}^{\mathrm{abs}},{\kappa }_{v,i}^{\mathrm{sca}}$ and g_{v, i} are adopted from the DSHARP opacity model (Birnstiel et al. 2018). The optical depth (838 μm) at the edge of photoevaporative cavity is ∼2.1 at t=0.60 Myr (see Fig. A.1). The τ_v values of the three bumps are approximately in the ratio of 4:2.4:1, whereas the corresponding dust surface density ratios are 3:2:1. The relatively high optical depth at the cavity edge arises not only from the larger dust surface density but also from the preferential removal of small particles by photoevaporation, which increases the fraction of particles larger than 1 mm. Since such large particles have higher opacities at the observation wavelength, the resulting τ_ν becomes significantly enhanced in this region.

Gárate et al. (2023) also studies the combined impact of early substructures and photoevaporation. They find that in disks with primordial substructures, once photoevaporation opens an inner cavity, further dust loss due to radial drift is largely prevented, and dust removal by wind entrainment has only a minor effect on the evolution of solids. Such disks can retain a substantial dust mass. The difference between their work and ours is that they do not include the planetesimal formation. Furthermore, they indicate that the ratio of the gap location relative to the disk’s characteristic radius r/R_c controls how much dust can be retained.

Fig. 6 Synthetic dust continuum images at 880 μm for three evolutionary stages of the protoplanetary disk model with three local pressure bumps at 30 au and 60 au in Sect. 4.3. Panel a: t=0.38 Myr, pre-cavity stage; Panel b: t=0.60 Myr, post-cavity stage when the photoevaporative cavity moves to 18 au; Panel c: t=0.80 Myr, when the photoevaporative cavity merges with the inner bump at 30 au. The blanks within the inner disks in Panels b) and c) correspond to the photoevaporative cavities.

4.4 Caveats

We point out several limitations in our study. First, we neglect the interaction of dust onto gas, commonly referred to as dust-to-gas backreaction, in all of our simulations. This effect can be pronounced when the solid density approaches that of the surrounding gas (ε=ρ_d/ρ_g ∼1). From the classical SI perspective, the drag force exerted by dust particles accelerates the gas’s azimuthal motion. This, in turn, suppresses the inward drift of solids and further increases the local dust-to-gas ratio, making planetesimal formation more favorable.

Nevertheless, near the pressure maximum of the photoevaporation cavity, this backreaction plays a dual role: whereas it slows down the drift of dust, the gas simultaneously gains angular momentum, which acts to widen the pressure bump. This broadening would, on the other hand, limit the further dust concentration. Consequently, the overall impact of backreaction on the SI and planetesimal formation can be more complicated. We also conducted an additional simulation based on run_fid that included the dust-to-gas backreaction, implemented through the DustPyLib interface (Stammler et al. 2023). The result is overall similar in terms of final planetesimal mass, which decreases from 31.4 M_⊕ to 30.7 M_⊕.

The other limitation is not accounting for the effect of the ice line. Our model employs a static disk temperature profile and therefore does not account for the movement of the ice lines. Furthermore, we do not consider composition-dependent sublimation as pebbles drift across the ice lines of various species. The process of sublimation, followed by diffusion and recondensation, could locally enhance the solid abundance near the ice line location (Ida & Guillot 2016; Drążkowska & Alibert 2017; Schoonenberg & Ormel 2017). Consequently, the omission of the above mechanisms likely leads to an underestimation of planetesimal formation near ice lines. Future models should incorporate an evolving disk thermal structure and track the compositional evolution of pebbles to properly model the influence of ice lines.

Here we also do not consider the bouncing barrier. Laboratory and numerical studies indicate that the bouncing barrier typically prevents further coagulation once the dust aggregates reach millimeter-centimeter size (Zsom et al. 2010). Neglecting this effect may therefore overestimate the mass fraction of pebbles that can grow to the sizes relevant for triggering the SI.

Finally, we note that turbulence plays a more complex role in realistic disk environments than that assumed in Lim et al. (2024). Large-scale structures developed from magnetorotational or hydrodynamic instabilities can create local pressure maxima that efficiently trap solids and promote planetesimal formation. In particular, Eriksson et al. (2025) finds that vortices can efficiently trap drifting dust at their pressure maxima, leading to locally enhanced dust-to-gas ratios. The weak turbulence and low relative velocities within the central regions of vortices facilitate particle growth, and once planetesimals form, they can continue to grow through efficient pebble accretion. Likewise, Carrera et al. (2025) identifies a positive feedback between dust coagulation and the SI, where enhanced dust concentration reduces turbulence and drift, promoting further particle growth and a stronger SI. Moreover, Huang & Bai (2025) demonstrates that the interplay among the VSI and the SI can produce strong dust clumping within vortices in the outer disk regions even at moderately low dust-to-gas ratios. These effects are not considered in our current work.

5 Conclusions

In this work, we employed the DustPy code to study the dust coagulation, fragmentation, and radial drift, in a viscously accreting protoplanetary disk subject to stellar X-ray photoevaporation. Photoevaporation opens an inner cavity at orbital radii of a few au, causing the disk’s inner edge to expand radially outward as gas disperses. This process efficiently removes gas via photoevaporative winds and therefore largely enhances the local dust-to-gas ratio at the pressure bump near the cavity’s edge, which serves as a natural site to trigger the SI. We demonstrate that planetesimal formation via the SI, driven by the stellar X-ray photoevaporation, is robust during the disk dispersal phase. In our fiducial model, in total 31.4 M_⊕ of planetesimals form within the 100 au disk region. Furthermore, the exact choice of the SI criteria has little impact on the results.

We explore how the key physical processes regulate both the onset of the SI and the final masses of planetesimals. We find that the efficiency of SI-induced planetesimal formation is strongly regulated by the available dust reservoir and its retention during disk evolution. Higher metallicity retains dust mass more effectively, directly enhancing planetesimal formation. Higher viscosity induces faster gas replenishment from the outer disk region, delays cavity opening, reduces the overall dust reservoir, and thus hinders planetesimal formation. Higher stellar X-ray luminosity promotes earlier photoevaporative cavity opening while preserving dust at the cavity edge, thereby favoring planetesimal formation. Larger disks prolong dust retention, thereby extending the duration of planetesimal formation. By retaining dust more efficiently, substructured disks create favorable conditions for planetesimal formation. To summarize, our models suggest that planetesimal formation is favored in disks with larger sizes, substructures, higher metallicities, lower viscosities, higher dust fragmentation threshold velocities, and/or around stars with higher X-ray luminosities. Future work will focus on the growth and migration of these forming planetesimals (Fang et al. 2023; Lau et al. 2022, 2024, 2025) starting from a realistic birth mass distribution (Schäfer et al. 2017; Liu et al. 2020).

Acknowledgements

We thank the anonymous referee for constructive feedback on the manuscript. This work is supported by the National Key R&D Program of China (2024YFA1611803). B.L. acknowledges support from the startup grant of the Bairen program from Zhejiang University. HY is supported by the National Natural Science Foundation of China (Grant No. 12473067). J.D. is funded by the European Union under the European Union’s Horizon Europe Research & Innovation Program 101040037 (PLANETOIDS).

References

Appendix A Optical depth

We compute the optical depth τ_ν following: $${\tau }_v={\sum }_i({\kappa }_{v,i}^{\mathrm{abs}}+{\kappa }_{v,i}^{\mathrm{sca},\mathrm{eff}}){\Sigma }_{\mathrm{d}}^i,$$(A.1) where ${\kappa }_{v,i}^{\text{abs}}$ and ${\kappa }_{v,i}^{\text{sca,eff}}=(1-{\mathrm{g}}_{v,i}){\kappa }_{v,i}^{\text{sca}}$ are the absorption and effective scattering opacities of dust species i, and ${\Sigma }_{\mathrm{d}}^i$ is the dust surface density of dust species i. The values of ${\kappa }_{v,i}^{\mathrm{abs}},{\kappa }_{v,i}^{\mathrm{sca}}$ and g_v,i are adopted from the DSHARP opacity model (Birnstiel et al. 2018).

Fig. A.1 Radial profiles of the optical depth τv at 838 μm for the three pressure bumps at t=0.60 Myr.

The optical depth (838 μm) at the edge of photoevaporative cavity is ∼2.1 at t=0.60 Myr. The τ_v values of the other two bumps are approximately ∼1.2 at ∼34 au and ∼0.5 at ∼70 au. The optical depth of the outer bump at ∼70 au is consistent with the results reported by Dullemond et al. (2018), who shows that observed outer disk rings typically have optical depths of order ∼0.5. Stammler et al. (2019) attributes these moderate optical depths to ongoing planetesimal formation, in agreement with our findings. In contrast, our results indicate significantly higher optical depths at the photoevaporative cavity edge and in the inner bump.

Appendix B Dust entrainment

In this work, we adopt the dust entrainment model of Sellek et al. (2020) as our fiducial prescription for dust removal by photoevaporation. However, subsequent studies by Booth & Clarke (2021) and Burn et al. (2022) suggest that the entrained dust size may be significantly smaller, motivating an alternative dust entrainment model driven by internal photoevaporation: $$a_{\mathrm{ent},\text{int}}=\sqrt{\frac{8}{\pi }}\frac{{\dot{\Sigma }}_{\mathrm{PE},\mathrm{g}}}{{\rho }_{•}{\Omega }_{\mathrm{K}}}\frac{H_{\text{base}}}{z_{\text{base}}}{(1+\frac{z_{\text{base}}^2}{r^2})}^{3/2},$$(B.1) where H_base ≃H_g denotes the gas scale height at the wind base, and z_base is the height at which the photoevaporative wind is launched. We adopt z_base/H_base=4 as suggested by Booth & Clarke (2021). The maximum entrained size a_ent,int is limited to only a few microns.

To assess the impact of these different assumptions, we performed a set of comparative simulations using three models: the fiducial model (FID), the pure gas photoevaporation model (PG), and the dust entrainment model driven by internal photoevaporation (IDE).

Fig. B.1 a compares the final planetesimal mass obtained at the end of the simulations, whereas Fig. B.1b shows the time evolution of the total planetesimal mass M_plt for the three models.

Fig. B.1 Panel a) Total planetesimal mass Mplt for three models FID (purple), PG (orange), and IDE (cyan). Panel b) Time evolution of the planetesimal mass Mplt for three models.

We find that there is no significant difference between the PG and IDE models, either in the time evolution of M_plt or in the final planetesimal mass. This similarity arises because, in the IDE model, only dust particles with sizes of a few microns are efficiently entrained by the photoevaporative winds, leaving the bulk of dust mass largely unaffected. In contrast, the FID model, which allows the removal of larger dust particles, systematically underestimates the final planetesimal mass. We will improve the treatment of dust entrainment in future work.

All Tables

All Figures

Fig. 1 Time and radial distance evolution of different quantities in the fiducial model. Panel a: gas (dashed lines) and dust (solid lines) surface density. Panel b: Stokes number of the maximum-size dust particles (solid lines) and the mean Stokes number of the dust particles (dashed lines) at each radial bin. Maximum-size refers to the size of the particle whose surface density is at the peak of the mass bins hereafter. Panel c: midplane density ratio ε (solid lines) and metallicity Z (dashed lines).

In the text

	Fig. 2 Synthetic 880 μm continuum images from the fiducial model. Panel a: t=0.25 Myr, pre-cavity stage; Panel b: t=0.35 Myr, onset of photoevaporative cavity opening; Panel c: t=0.45 Myr, post-cavity stage.
In the text

Fig. 3 Evolution of dust mass (black lines), planetesimal mass (red lines), and dust mass fraction within the pressure bump in three Stokes number bins (light blue for St<10−2, blue for 10−2 ≤ St ≤ 10−1, and deep blue for St>10−1), across four evolutionary stages: (i) pre-cavity, (ii) cavity at 10 au, (iii) cavity at 20 au, and (iv) cavity at 100 au. These panels show how initial conditions affect dust retention and planetesimal formation: Panel a: fiducial model (Z=0.01, LX=2.6 LX, ⊙, Mdisk=0.05 M⊙, Rc=60 au, α=10−4) serves as the baseline for comparison, showing balanced dust evolution and planetesimal production across the disk. Panel b: higher metallicity (Z=0.02) facilitates the retention of dust mass, thereby promoting planetesimal formation. Panel c: decreased X-ray luminosity (from 2.6 LX, ⊙ to 1.0 LX, ⊙) diminishes dust mass and consequently impedes planetesimal formation. Panel d: higher viscosity parameters (α=10−3, δt=10−4) induce gas replenishment from the outer disk, delay cavity opening, reduce the overall dust reservoir, and thus hinder planetesimal formation. Panel e: a more massive and extended disk (Mdisk=0.1 M⊙, Rc=120 au) prolongs dust retention, thereby extending the duration of planetesimal formation. Panel f: a less massive but extended disk (Mdisk=0.05 M⊙, Rc=120 au) lowers the inner-disk dust density, causing photoevaporation to open a cavity earlier. This makes more dust retained at the onset of photoevaporative cavity, which also provides more building blocks for planetesimal formation.

In the text

	Fig. 4 Distribution of dust as a function of radial distance and particle size at t=0.30 Myr for run_fid (left) and run_frag (right). The color bar represents the surface density of dust per logarithmic mass interval. The green line is the drift limit (Eq. (14)), whereas the blue line is the fragmentation limit (Eq. (15)). The four white lines correspond to St=[10−3, 10−2, 10−1, 100].
In the text

	Fig. 5 Comparison of the radial mass distributions of planetesimals at different disk regions. The curves show the results for models with varying parameters: run_fid (black), run_metal (blue), run_lumi (orange), run_alpha (purple), run_dsize1 (red), and run_dsize2 (cyan).
In the text

Fig. 6 Synthetic dust continuum images at 880 μm for three evolutionary stages of the protoplanetary disk model with three local pressure bumps at 30 au and 60 au in Sect. 4.3. Panel a: t=0.38 Myr, pre-cavity stage; Panel b: t=0.60 Myr, post-cavity stage when the photoevaporative cavity moves to 18 au; Panel c: t=0.80 Myr, when the photoevaporative cavity merges with the inner bump at 30 au. The blanks within the inner disks in Panels b) and c) correspond to the photoevaporative cavities.

In the text

	Fig. A.1 Radial profiles of the optical depth τv at 838 μm for the three pressure bumps at t=0.60 Myr.
In the text

	Fig. B.1 Panel a) Total planetesimal mass Mplt for three models FID (purple), PG (orange), and IDE (cyan). Panel b) Time evolution of the planetesimal mass Mplt for three models.
In the text