© The Authors 2025

1 Introduction

The formation of planetesimals, typically hundreds of kilometers in diameter, from micrometer-sized dust particles is a fundamental stage in planetary evolution according to the core accretion theory (Pollack et al. 1996). A significant challenge in this process is generally known as the “meter-size barrier”, which refers to the difficulty in growing dust particles beyond one meter due to rapid radial drift and destructive collisions (Blum & Wurm 2008). One potential solution to overcome this barrier is the presence of dust traps – regions where the aerodynamic drag force on dust grains weakens or vanishes entirely, facilitating dust accumulation and growth. Pressure bumps, a typical kind of dust trap arising from a local maximum in gas pressure, have received increasing attention in recent years as a favorable site for planet formation (e.g., Morbidelli 2020; Guilera et al. 2020; Chambers 2021; Jiang & Ormel 2023; Sándor et al. 2024; Lau et al. 2022, 2024). Since dust particles drift toward higher pressure (Weidenschilling 1977; Nakagawa et al. 1986; Takeuchi & Lin 2002), pressure bumps act as barriers to inward drifting dust, accumulating the dust from the outer disk and increasing dust concentration within the bump. The locally enriched dust-to-gas ratio could trigger streaming instability (Youdin & Goodman 2005; Youdin & Lithwick 2007; Johansen et al. 2007, 2009; Bai & Stone 2010), leading to the formation of planetesimals. The compositional and dynamical properties of planetesimals that formed through streaming instability are generally consistent with those of comets and Kuiper Belt objects (Blum et al. 2017; Nesvorný et al. 2019).

Recent high-resolution interferometry observations from the Atacama Large Millimeter/submillimeter Array (ALMA) have shown that substructures are prevalent in protoplanetary disks. These substructures primarily appear as axisymmetric rings, as documented in various surveys (e.g., Andrews et al. 2018; Long et al. 2018; Cieza et al. 2021). Although these observations focus on the largest and brightest disks, comparisons between observational data of disk populations and theoretical models imply that such substructures must be common even in unresolved disks (Toci et al. 2021; Zormpas et al. 2022; Delussu et al. 2024). A detailed study by Dullemond et al. (2018) demonstrates that dust trapping in local pressure maxima is consistent with the rings sampled from the DSHARP survey. Kinematic studies with ALMA by Teague et al. (2018a,b) further confirm that dust rings are associated with pressure bumps in the analyzed disks. Stammler et al. (2019) propose that planetesimal formation through streaming instability, where the midplane dust-to-gas ratio is regulated, could account for the observed similarities in the optical depths across the DSHARP rings identified by Dullemond et al. (2018). Hydrodynamical simulations with self-gravity by Carrera et al. (2021) indicate that even a small pressure bump could initiate planetesimal formation via streaming instability with centimeter-sized grains. However, this process may be problematic for millimeter-sized particles (Carrera & Simon 2022), which are most accessible to (sub)millimeter observations with ALMA.

The origin of pressure bumps, however, remains uncertain. The most widely proposed mechanism is planet-disk interaction, where a massive planet opens a gap in the disk, resulting in a pressure bump just outside the gap (Rice et al. 2006; Pinilla et al. 2012; Dipierro et al. 2015; Dong et al. 2017; Zhang et al. 2018). However, among the many disks observed by direct imaging, PDS 70 remains the only system where accreting exoplanets have been unambiguously identified (Keppler et al. 2018; Haffert et al. 2019). Furthermore, the formation of such massive planets themselves needs to be explained independently of planet- induced pressure bumps. Therefore, non-planetary mechanisms that might induce pressure bumps should also be considered, including dead zone boundaries (Flock et al. 2015; Pinilla et al. 2016), zonal flows (Bai & Stone 2014), snowlines (Kretke & Lin 2007; Saito & Sirono 2011), photoevaporation (Owen & Kollmeier 2019), and infall of material (Kuznetsova et al. 2022).

Recent discoveries of infall streamers indicate that infall can occur at various evolutionary stages of protostars, including Class 0 protostars, which are embedded in dense gas and dust envelopes (e.g., Tobin et al. 2012; Tokuda et al. 2018; Pineda et al. 2020; Thieme et al. 2022; Murillo et al. 2022); Class I protostars, where protoplanetary disks are gradually forming (e.g., Yen et al. 2014; Garufi et al. 2022; Valdivia-Mena et al. 2022; Flores et al. 2023; Cacciapuoti et al. 2024); and Class II protostars, where protoplanetary disks are fully developed (e.g., Tang et al. 2012; Alves et al. 2020; Ginski et al. 2021; Huang et al. 2021; Gupta et al. 2023). Winter et al. (2024) estimate that 20% to 70% of protoplanetary disks accrete more than 50% of their material through late-stage infall. While numerous works have modeled the formation and evolution of planetary cores in pressure bumps caused by various mechanisms (e.g., Morbidelli 2020; Guilera et al. 2020; Chambers 2021; Jiang & Ormel 2023; Sándor et al. 2024; Lau et al. 2022, 2024), as of yet no study has simulated planet formation specifically in infall-induced pressure bumps.

In this work, we present a numerical model for the formation of planetesimals in an axisymmetric pressure bump induced by late-stage infall onto a protoplanetary disk. Section 2 describes the physical models for the structure and evolution of gas and dust in the disk, the infall process, and the formation of planetesimals from the accumulated dust. Section 3 outlines our simulation setup using the gas and dust evolution code DustPy (Stammler & Birnstiel 2022), along with the parameter space explored. The simulation results and analysis are detailed in Section 4. We provide discussions of our model and results in Section 5, and conclude the study in Section 6.

2 Model

Our model consists of disk evolution, infall, and planetesimal formation. We assume a 1D axisymmetric protoplanetary disk orbiting around a Solar-mass central star, where the gas undergoes viscous evolution and the dust drifts and grows.

2.1 Disk structure and evolution

2.1.1 Gas component

The initial gas surface density profile of the protoplan- etary disk adopts the self-similar solution derived by Lynden-Bell & Pringle (1974), $${\text{Σ}}_{\text{g,init\hspace{0.17em}}}=\frac{M_{\text{disk\hspace{0.17em}}}}{2\pi r_{\text{c}}^2}{\left(\frac{r}{r_{\text{c}}}\right)}^{-1}\exp \left(-\frac{r}{r_{\text{c}}}\right),$$(1)

where r is the distance from the central star, M_disk is the initial mass of the disk, and r_c is the characteristic radius of the disk. Assuming M_disk = 0.05 M_⊙ and r_c = 100 au, for instance, it yields Σ_g,init ≈ 700 g cm⁻² at r = 1 au.

The midplane gas temperature is calculated assuming a passively irradiated disk with a constant irradiation angle of 0.05 (Chiang & Goldreich 1997; Dullemond et al. 2001), $$T_{\text{g}}\approx 221.3{\left(\frac{r}{\text{au}}\right)}^{-1/2}\text{\hspace{0.17em}K},$$(2)

from which the isothermal sound speed at the midplane $c_{\text{s}}=\sqrt{k_{\text{B}}{T}_{\text{g}}/\mu m_{\text{p}}}$ can be calculated with the Boltzmann constant k_B, the mean molecular weight of the gas µ = 2.3, and the proton mass m_p.

The protoplanetary disk is assumed to be in vertical hydrostatic equilibrium. Given the pressure scale height of the gas defined as H_g ≡ c_s/Ω_K, where Ω_K is the Keplerian orbital frequency, the aspect ratio of the disk is quantified by $$h=\frac{H_{\text{g}}}{r}\approx 0.03{\left(\frac{r}{\text{au}}\right)}^{1/4}.$$(3)

The midplane gas volume density can be derived according to vertical hydrostatic equilibrium as $${\rho }_{\text{g}}=\frac{{\text{Σ}}_{\text{g}}}{\sqrt{2\pi }{H}_{\text{g}}}.$$(4)

The midplane gas pressure is then given by the ideal gas law $P={\rho }_{\text{g}}{c}_{\text{s}}^2$, from which we can obtain the midplane pressure gradient parameter $$\eta =-\frac{h^2}{2}\frac{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}P}{\partial \text{\hspace{0.17em}ln\hspace{0.17em}}r}.$$(5)

This quantity characterizes the degree of “sub-Keplerianity” of the gas, given by the relation $v_{\text{g},\phi }^2=(1-2\eta )v_{\text{K}}^2$, where v_{g,φe flux follows a Gaussian distribution} is the azimuthal velocity of the gas and v_K is the Keplerian orbital velocity. For small values of η, this expression simplifies to v_g,φ = (1 − η)v_K.

The evolution of the gas follows the viscous advectiondiffusion equation (Lüst 1952; Lynden-Bell & Pringle 1974; Pringle 1981), $$\frac{\partial {\text{Σ}}_{\text{g}}}{\partial t}=\frac{3}{r}\frac{\partial }{\partial r}\left[\sqrt{r}\frac{\partial }{\partial r}\left(\upsilon {\text{Σ}}_{\text{g}}\sqrt{r}\right)\right]+S_{\text{g},\text{\hspace{0.17em}infall\hspace{0.17em}}},$$(6)

with the kinematic viscosity ν = αc_s H_g adopting the α prescription introduced by Shakura & Sunyaev (1973), and the external source term S_g,infall referring to the flux of the infalling gas. Neglecting the dynamic back reaction of dust onto the gas, the gas radial velocity due to viscous accretion is given by $${\upsilon }_{\text{g},r}=-\frac{3}{{\text{Σ}}_{\text{g}}\sqrt{r}}\frac{\partial }{\partial r}\left(\upsilon {\text{Σ}}_{\text{g}}\sqrt{r}\right).$$(7)

2.1.2 Dust component

The initial surface density of the dust is assumed to be proportional to that of the gas, $${\text{Σ}}_{\text{d},\text{\hspace{0.17em}init\hspace{0.17em}}}={ϵ}_{\text{init\hspace{0.17em}}}{\text{Σ}}_{\text{g},\text{\hspace{0.17em}init\hspace{0.17em}}},$$(8)

where the initial global dust-to-gas ratio ϵ_init is set to the Solar metallicity of 0.01. The initial size distribution of the dust particles n(a) follows the MRN (Mathis-Rumpl-Nordsieck) size distribution of interstellar grains introduced by Mathis et al. (1977), expressed as n(a) ∝ a^−7/2 up to the maximum initial particle radius of 1 µm. A bulk mass density of ρ_s = 1.67 g cm⁻³ is assumed. The Stokes number St_i of each particle mass bin m_i (see Section 3 for details on particle mass binning) is determined by considering both the Epstein and Stokes I drag regimes.

The dust scale height of each mass bin H_d,i is derived by solving a vertical settling-mixing equilibrium equation by Dubrulle et al. (1995), calculated as $$H_{\text{d},i}=H_{\text{g}}\sqrt{\frac{{\delta }_{\text{vert}}}{{\delta }_{\text{vert}}+{\text{St}}_i}},$$(9)

where the vertical mixing parameter δ_vert = α is assumed. Then, the midplane dust volume density of each mass bin ρ_d,i is obtained via vertical hydrostatic equilibrium.

The evolution of the dust surface density Σ_d,i involves dust transport, dust growth, and external sources from infall. Dust transport is described by the advection-diffusion equation (Clarke & Pringle 1988; Birnstiel et al. 2010), $$\frac{\partial {\text{Σ}}_{\text{d},i}}{\partial t}+\frac{1}{r}\frac{\partial }{\partial r}\left[r{\text{Σ}}_{\text{d},i}{v}_{\text{d},i}-r{D}_{\text{d},i}{\text{Σ}}_{\text{g}}\frac{\partial }{\partial r}\left(\frac{{\text{Σ}}_{\text{d},i}}{{\text{Σ}}_{\text{g}}}\right)\right]=S_{\text{d},\text{infall},i}$$(10)

where $D_{\text{d},i}={\delta }_{\text{rad}}{c}_{\text{s}}{H}_{\text{g}}/\left(1+{\text{St}}_i^2\right)$ is the dust diffusivity (Youdin & Lithwick 2007) assuming the radial mixing parameter δ_rad = α, and the external source term S_d,infall,i is the dust infall flux. The radial velocity of the dust v_d,i, according to Nakagawa et al. (1986) and Takeuchi & Lin (2002), is given by $${\upsilon }_{\text{d},i}\simeq -\frac{2\eta {\upsilon }_{\text{K}}}{{\text{St}}_i+{\text{St}}_i^{-1}}.$$(11)

Dust growth by coagulation is calculated by solving the Smolu- chowski equation (Smoluchowski 1916).

The radial drift speed of the dust rises as the particle sizes increase within St < 1, thus setting a limit for dust particle sizes as the drift rate exceeds the coagulation rate. Another limit is fragmentation, which occurs when the impact velocity between two dust grains exceeds a critical fragmentation velocity v_frag. Laboratory experiments have measured v_frag around 100 cm s⁻¹ for silicate grains (Blum & Wurm 2008; Güttler et al. 2010). Water-ice particles have a significantly higher v_frag of over 1000 cm s⁻¹ due to their approximately tenfold higher surface energies, as confirmed experimentally by Gundlach & Blum (2015). The disk population study by Delussu et al. (2024) suggests a range of v_frag ≥ 500 cm s⁻¹. In this work, we adopt a relatively conservative range, limiting v_frag to values below 500 cm s⁻¹. For further details on the gas and dust evolution algorithm, readers can refer to Stammler & Birnstiel (2022).

Fig. 1 Illustration of the infall model, showing profiles of gas surface density (blue curves) and pressure (orange curves). The dashed curves represent the initial profiles, while the solid curves depict the profiles after infall. The horizontal blue and orange arrows indicate the directions of gas and dust radial velocities, respectively.

2.2 Infall model

Gas and dust infall onto the protoplanetary disk over a duration of t_infall, starting at t = 0. The infall rate is constant during this process. We define the infall rate of gas by Ṁ_infall ≡ M_infall /t_infall, where M_infall is the total mass of the infalling gas. The global dust-to-gas ratio of the infalling materials equals the initial condition of the disk, ϵ_init = 0.01, thus the total mass of the infalling materials consisting of gas and dust is 1.01 M_infall. The initial particle sizes of the infalling dust also follow the MRN size distribution.

The infalling materials are distributed axisymmetrically in the disk, and the flux follows a Gaussian distribution over r, $$S_{\text{g},\text{\hspace{0.17em}infall\hspace{0.17em}}}=\frac{{\dot{M}}_{\text{infall\hspace{0.17em}}}}{A}\exp \left[-\frac{{\left(r-r_{\text{infall\hspace{0.17em}}}\right)}^2}{2{\sigma }_r^2}\right].$$(12)

where r_infall is the central radius of infall, assumed to be fixed over time, and σr describes the width of the Gaussian infall. The normalizing area factor A ensures the integral of S_g,infall over the disk area yields Ṁ_infall, formulated as $$A={\displaystyle {\int }_{r_{\min }}^{r_{\max }}2}\pi r\exp \left[-\frac{{\left(r-r_{\text{infall\hspace{0.17em}}}\right)}^2}{2{\sigma }_r^2}\right]\text{d}r,$$(13)

with the inner and outer boundaries of the disk r_min and r_max.

Figure 1 provides an illustration of the infall model, with the central infall radius at r_infall = 40 au as indicated by the black arrow. The infall width parameter is set to σ_r = 4 au, and gas with a total mass of 10⁴ M_⊕ falls onto a disk with an initial mass of M_disk = 0.05 M_⊙ and characteristic radius r_c = 100 au. A local maximum in gas surface density, marked by the blue dot, appears slightly inward of r_infall. Since the sound speed c_s decreases with radial distance, a local maximum in pressure forms even closer to the star, as indicated by the orange dot. This pressure bump results in a reversal of the pressure gradient across it, distinguishing sub-Keplerian (green-shaded region, η > 0) and super-Keplerian (red-shaded region, η < 0) gas. The orange horizontal arrows illustrate the dust drift direction, which follows the pressure gradient as given by Eq. (11). Dust particles are drawn toward the pressure peak, leading to increased dust concentration at this location. The gas radial velocity, calculated from Eq. (7), changes direction at the vertical dotted blue line. Over time, viscous diffusion gradually smooths the pressure bump, making the dust trap a temporary feature. Because this boundary lies beyond the local maxima in gas surface density and pressure, viscous spreading causes the pressure maximum to shift inward. Since this work neglects the back reaction of dust on gas, which generally causes the gas to flow toward lower pressure (Nakagawa et al. 1986), the pressure gradient does not contribute to smoothing the bump.

2.3 Planetesimal formation

As the dust accumulates and grows at the pressure bump induced by infall, the midplane dust-to-gas ratio gradually increases, and streaming instability can be triggered to transform dust into planetesimals. The prescriptions for the transformation are based on the methods in Drążkowska et al. (2016) and Schoonenberg et al. (2018), which considered a midplane dust-to-gas ratio, $${ϵ}_{\text{mid}}=\frac{{\rho }_{\text{d}}}{{\rho }_{\text{g}}}={\displaystyle \sum _i\frac{{\rho }_{\text{d},i}}{{\rho }_{\text{g}}}},$$(14)

exceeding unity as the condition for streaming instability (Drążkowska & Dullemond 2014). In addition, following the prescriptions in Lau et al. (2024), the condition for gravitational collapse of a dense filament induced by streaming instability is determined by a Toomre-like criterion, requiring Q_p < 1. The parameter Q_p is given by $$Q_{\text{p}}=\sqrt{\frac{\delta }{{\text{St}}_{\text{avg\hspace{0.17em}}}}}\frac{c_{\text{s}}{\mathrm{\Omega }}_{\text{K}}}{\pi G{\text{Σ}}_{\text{d},\text{\hspace{0.17em}local\hspace{0.17em}}}}=\sqrt{\frac{\delta }{{\text{St}}_{\text{avg}}}\frac{Q}{{\text{Σ}}_{\text{d},\text{\hspace{0.17em}local\hspace{0.17em}}}/{\text{Σ}}_{\text{g}}}}$$(15)

as defined in Gerbig & Li (2023). Here, Q = c_sΩ_K /(πGΣ_g) is the standard Toomre parameter (Toomre 1964). The term δ represents the small-scale diffusion parameter. The density-averaged Stokes number is defined as ${\text{St}}_{\text{avg}}={\displaystyle \sum _i{\text{St}}_i}{\text{Σ}}_{\text{d},i}/{\displaystyle \sum _i{\text{Σ}}_{\text{d},i}}$. The local dust surface density at the dense filament, Σ_d,local, is assumed to be ten times Σ_d. According to the measurements in streaming instability simulations by Schreiber & Klahr (2018), the value of δ ranges from 10⁻⁶ to 10⁻⁴. Since a lower δ is beneficial to planetesimal formation, we fixed δ at a conservative value of 10⁻⁴ to make planetesimal formation more challenging.

Once both conditions, ϵ_mid ≥ 1 and Q_p< 1, are satisfied simultaneously in a radial cell (see Section 3 for details on the radial grid), the dust surface density in that cell is transformed into the surface density of planetesimals by $$\frac{\partial {\text{Σ}}_{\text{plts}}}{\partial t}=-{\displaystyle \sum _i{\left(\frac{\partial {\text{Σ}}_{\text{d},i}}{\partial t}\right)}_{\text{plts}}},$$(16)

with the reduction rate of dust surface density of each dust species i that contributes to planetesimal masses $${\left(\frac{\partial {\text{Σ}}_{\text{d},i}}{\partial t}\right)}_{\text{plts}}=-{\mathcal{P}}_{\text{pf}}\frac{\zeta {\text{Σ}}_{\text{d},i}}{t_{\text{sett},i}}=-{\mathcal{P}}_{\text{pf}}\zeta {\text{Σ}}_{\text{d},i}{\text{St}}_i{\mathrm{\Omega }}_{\text{K}},$$(17)

The settling time t_sett,i = (St_i Ω_K)⁻¹ is the timescale on which streaming instability filaments form (Yang et al. 2017). The planetesimal formation efficiency ζ describes the fraction of dust transformed into planetesimals per settling time, set at 10⁻³. The parameter ζ is not yet well understood; however, it does not impact the conditions required to trigger planetesimal formation. The term 𝒫_pf is a smooth activation function designed to prevent abrupt initiation of planetesimal formation in Miller et al. (2021), where it was defined as a hyperbolic tangent function of ϵ_mid. Unlike Miller et al. (2021), here it is defined as a function of Q_p, formulated as $${\mathcal{P}}_{\text{pf}}=\frac{1}{2}\left[1-\tanh \left(\frac{Q_{\text{p}}-0.75}{n}\right)\right],$$(18)

with the smoothness parameter n. We set n = 0.2 and made the function center around Q_p = 0.75 so that the value at the planetesimal formation threshold Q_p = 1 is positive but rather low, as shown in Figure 2. Each of our simulations eventually produces a surface density distribution of planetesimals.

In the outer regions of the disk, the relative velocity between dust grains can remain consistently below vfrag, preventing fragmentation from occurring. If the pressure bump is located in this region, some dust particles will grow to St >> 1 in the absence of a fragmentation limit. However, streaming instability requires dust grains that are moderately coupled to the gas (Youdin & Goodman 2005; Johansen et al. 2007; Bai & Stone 2010), meaning that excessively large dust particles will not participate in streaming instability. Numerous studies have shown that when the vertically integrated dust-to-gas ratio ϵ is around the Solar metallicity (∼0.01), strong clumping induced by streaming instability can only occur for dust grains with 0.01 ≲ St ≲ 1 (Carrera et al. 2015; Yang et al. 2017; Li & Youdin 2021; Simon et al. 2022; Birnstiel 2024). Therefore, all calculations from Eqs. (15) to (18) only consider dust species with St < 1. A lower limit for St is unnecessary because the mass fraction of dust species with St ≪ 0.01 is negligible once the dust trap has formed.

Fig. 2 Smooth activation functions in terms of Qp with smoothness parameter n = 0.2, centered at Qp= 0.75. Gravitational collapse occurs when Qp < 1.

3 Simulation setup and parameter space

Our simulations were performed using the code DustPy¹ (Stammler & Birnstiel 2022). DustPy is a Python package to simulate the radial evolution of gas and dust in protoplanetary disks, including viscous gas evolution, dust advection and diffusion, as well as dust growth by coagulation and fragmentation, based on the model of Birnstiel et al. (2010).

The dust particle masses are binned logarithmically from 10⁻¹² to 10⁸ g, using a total of 141 bins. The radial grid spans from r_min = 1 au to r_max = 10³ au, with refinement around the infall region as follows: First, a logarithmic radial grid from r_min to r_max with 100 cells is created. Then, the section of the grid between (r_infall − 3σ_r) and (r_infall + 3σ_r) is replaced with a logarithmic grid consisting of 60 cells, where 99.7% of the infall materials are dumped.

There are eight parameters considered in our model: the gas viscosity parameter α, the dust fragmentation velocity v_frag, the initial disk mass M_disk, the characteristic disk radius r_c, the central infall radius r_infall, the gas infall rate Ṁ_infall, the infall duration t_infall, and the relative width of the infall region σ ≡ σ_r/r_infall. The value ranges of the parameters are listed in Table 1. Infall begins at the start of the simulations, and all simulations terminate at t_tot = 500 kyr, which is long enough compared with t_infall to determine whether planetesimals can form.

The range of α values aligns with those adopted in many recent simulation studies on pressure bumps (e.g., Chambers 2021; Lau et al. 2024; Sándor et al. 2024). Observational estimates (e.g., Dullemond et al. 2018; Doi & Kataoka 2021; Jiang et al. 2024) indicate that α ∼ 10⁻⁴ is consistent with the majority of observed disks. The value ranges of the disk parameters M_disk and r_c align with current observational measurements (Andrews 2020). The value ranges of the infall parameters Ṁ_infall and t_infall are determined according to the measurements in recent observations of young stellar objects with infall streamers (Alves et al. 2020; Pineda et al. 2020; Cacciapuoti et al. 2024). For instance, Cacciapuoti et al. (2024) observe a 4000 au streamer infalling onto the disk of a Class I protostar ([MGM2012] 512), whose infall rate is measured to be 1.5 × 10⁻⁶ M_⊙ yr⁻¹ (≈ 0.5 M_⊕ yr⁻¹) with a free-fall timescale of 50 kyr. In our case, assuming a streamer with a free-fall timescale of $t_{\text{ff}}=\sqrt{R^3/G{M}_{\star }}=50$ kyr and given M_⋆ = M_⊙, the length of the streamer is estimated to be R ≈ 4600 au, which is consistent with the order of magnitude observed.

4 Results

We present the detailed processes of pressure bump formation, dust accumulation, and planetesimal formation using a fiducial setup. The impact of each individual parameter on gas properties and planetesimal mass evolution is analyzed using the one- factor-at-a-time method. The effects of these parameters on the success of planetesimal formation are statistically investigated through random sampling.

4.1 Fiducial setup

The fiducial values of the parameters are displayed in Table 1. In the fiducial simulation, planetesimals begin to form at t = 1.12 × 10⁵ yr. Figure 3 shows the radial and particle mass distribution of the dust at this moment. The dust surface density σ_d here is defined to be independent of the particle mass grid (Stammler & Birnstiel 2022). The drift limit, denoted by the green contour, is reached when the dust grains grow to a size where the drift timescale becomes shorter than the growth timescale (Birnstiel et al. 2012), given by $${\text{St}}_{\text{drift}}=\frac{ϵ}{h^2}{\left|\frac{\partial \text{\hspace{0.17em}ln}P}{\partial \text{\hspace{0.17em}ln}r}\right|}^{-1},$$(19)

where ϵ = Σ_d/Σ_g is the vertically integrated dust-to-gas ratio. The two spikes on this curve indicate the positions where η changes sign. The fragmentation limit, denoted by the blue contour, is an approximation considering only turbulent motion, calculated as $${\text{St}}_{\text{frag}}=\frac{3-\sqrt{9-4b^2}}{2b},$$(20)

with $b=v_{\text{frag\hspace{0.17em}}}^2/\alpha c_{\text{s}}^2$. For b << 1 (i.e., St_frag << 1), this expression simplifies to St_frag = b/3. Dust growth is limited mainly by the fragmentation limit at the inner part of the disk, and by the drift limit at the outer part. The pressure gradient parameter η shown in Figure 4 has switched sign between ∼26 and ∼36 au by this time, forming a pressure bump that peaks at ∼36 au. The negative part of η forces the dust in this super-Keplerian region to drift outward, traps the infalling dust, and stops the dust from the outer disk drifting inward, causing strong dust accumulation at the peak of the pressure bump. The increasing dust-to-gas ratio in the pressure bump accelerates dust growth, because the dust growth timescale given by (Birnstiel et al. 2012) $$t_{\text{grow\hspace{0.17em}}}\simeq \frac{1}{ϵ{\mathrm{\Omega }}_{\text{K}}}$$(21)

is shortened. As shown in Figure 3, a significant proportion of the dust in the pressure bump grows to St > 0.1, whose growth is constrained by the fragmentation limit. The midplane dust-togas ratio ϵ_mid gradually rises and reaches the threshold 1, while the increasing St_avg and Σ_d,local reduce Q_p below the threshold 1, as Figure 4 shows. Then, streaming instability is triggered in the pressure bump, transforming the surface density of the dust into that of planetesimals.

Planetesimal formation stops at t = 2.74 × 10⁵ yr. Figure 5 shows the radial profiles of the gas, dust, and planetesimal surface densities as well as the midplane dust-to-gas ratio ϵ_mid, at the start and end of the planetesimal formation process, respectively. The infalling gas produces a local maximum slightly inside r_infall = 40 au in the gas surface density profile, which induces a pressure bump slightly inside the local maximum of the gas surface density Σ_g. By the time when planetesimals start to form, the dust accumulated by the pressure bump has formed a sharp peak in the dust surface density profile, creating a similar shape for the profile of ϵ_mid. At this moment, planetesimal formation is initiated in only one radial cell where ϵ_mid reaches 1, which is in reality a thin planetesimal ring in the protoplan- etary disk. During the period of planetesimal formation, the viscous evolution of the gas gradually smooths out the bump of Σ_g and makes it expand, and meanwhile the peaks of the pressure bump, the dust surface density Σ_d, and ϵ_mid are shifted inward. At t = 2.42 × 10⁵ yr, η becomes positive throughout the disk, implying that the pressure trap vanishes. Then 32 kyr later, the maximum of ϵ_mid drops below 1, thus the planetesimal formation process eventually ceases. At this point, a planetesimal belt with a width of ∼8 au has formed, spreading from ∼36 to ∼28 au as a result of the viscous diffusion of the pressure bump. This process closely resembles the simulation results of Miller et al. (2021), where a migrating planet creates a moving dust trap that forms a planetesimal belt.

Figure 6 shows the time evolution of the gas, dust, and planetesimal masses over the entire simulation. During the first 50 kyr, both the masses of the gas and dust rise due to the infall process. Over the remaining 450 kyr, the gas mass decreases slowly because of viscous accretion, while the loss of the dust mass is contributed by radial drift and transformation into planetesimals. The rapid decline of the dust mass over the last 100 kyr is because the massive dust trapped in the pressure bump that is not converted into planetesimals drifts inward after the pressure bump vanishes. Planetesimals with a total mass of 9.97 M_⊕ formed over 1.62 × 10⁵ years, which accounts for 3.7% of the total dust mass (the initial dust mass in the disk plus the infalling dust mass). A movie showing the result of the fiducial simulation is available online.

Fig. 3 Dust distribution of the fiducial simulation at the start of planetesimal formation. The dust surface density σd is independent of the mass grid. The green contour indicates the drift limit, while the blue contour marks the fragmentation limit. White contours represent particle masses corresponding to Stokes numbers St = 1 (solid) and St = {10−1, 10−2, 10−3 } (dotted), respectively.

Fig. 4 Midplane dust-to-gas ratio ϵmid, Toomre-like value Qp, and pressure gradient parameter η around the infall region of the fiducial simulation at the start of planetesimal formation. The vertical dotted gray line indicates rinfall = 40 au. The horizontal dashed gray line represents the critical conditions for streaming instability, where ϵmid = 1 and Qp = 1, as well as the condition for the presence of a pressure bump, η = 0.

Fig. 5 Radial surface density distributions of gas, dust, and planetesimals, as well as the midplane dust-to-gas ratio ϵmid for the fiducial simulation at the start (top) and end (bottom) of planetesimal formation. The vertical dotted gray line indicates rinfall = 40 au. The horizontal dashed gray lines represent ϵmid = 1.

Fig. 6 Mass evolution of the gas, dust, and planetesimals over 500 kyr in the fiducial simulation. Dotted lines indicate the start and end of planetesimal formation.

Fig. 7 Gas surface density distributions at t = 100 kyr for varying parameters compared with the fiducial setup. Each panel is labeled with the varied parameter, with green and orange curves representing decreased and increased parameter values in Table 1, respectively.

Fig. 8 Pressure gradient parameter η at t = 100 kyr for varying parameters compared with the fiducial setup. Each panel is labeled with the varied parameter, with green and orange curves representing decreased and increased parameter values in Table 1, respectively. The horizontal dashed gray lines represent η = 0.

4.2 One-factor-at-a-time analysis

In each simulation, only one parameter – except for v_frag (as it does not affect gas properties) – is varied to its maximum or minimum value from the ranges listed in Table 1, while all other parameters are held at their fiducial values. The resulting distributions of gas surface density and the midplane pressure gradient parameter η at 100 kyr for these simulations are illustrated in Figures 7 and 8, respectively. Figure 9 presents the evolution of planetesimal mass in the simulations where planetesimal formation occurs. We select t = 100 kyr as the snapshot time because, by this point, the infall stages have concluded across all simulations. Furthermore, in most simulations where planetesimals form (except when r_infall = 80, au), the trigger time for planetesimal formation occurs around t = 100 kyr, as shown in Figure 9. In simulations that do not produce planetesimals, this time is close to the point when their ϵ_mid values reach their maximum, though still below 1. Table 2 summarizes the simulation outcomes, including the start and end times of planetesimal formation, the total planetesimal masses formed, and the conversion efficiencies from dust to planetesimals.

For the viscosity parameter α, the results for α = 1 × 10⁻⁴ (green) and α = 4 × 10⁻⁴ (orange) are presented. With the same amount of gas dumped onto identical disks, the Gaussian-like gas surface density profile is broader and lower for higher α compared to lower α. The pressure gradient parameter η is correlated with the slope of the gas surface density profile by the relations in Eqs. (4) and (5). Consequently, the region with negative η spans a narrower range of r for higher α, and has lower absolute values of η. For α = 1 × 10⁻⁴, the planetesimal formation process is much longer-lasting and more productive than that of the fiducial setup. In contrast, no planetesimal is formed in the simulation with α = 4 × 10⁻⁴ since ϵ_mid never reaches 1, though the condition for Q_p can be satisfied.

For the initial disk mass M_disk, the results for M_disk = 0.02 M_⊙ (green) and M_disk = 0.1 M_⊙ (orange) are presented. Although the amount of infalling gas is the same, the Gaussian- like gas surface density profiles differ in height due to the varying initial profiles of the disks. For smaller M_disk, the bump in the gas surface density profile is more prominent, leading to a steeper slope and consequently a higher and broader absolute negative pressure gradient. For M_disk = 0.02 M_⊙, the planetesimal formation process is longer-lasting and more productive compared to the fiducial setup. In contrast, no planetesimal is formed for M_disk = 0.1 M_⊙ since ϵ_mid never reaches 1, though the condition for Q_p can be satisfied.

For the characteristic radius of disk r_c, the results for r_c = 50 au (green) and r_c = 200 au (orange) are presented. Since r_infall = 40 au is always smaller than r_c here, the initial gas surface density is higher in the infall region for smaller r_c, resulting in a flatter gas surface density bump. Hence, the region with negative η is broader and has higher absolute values of η for larger r_c. The planetesimal mass produced with r_c = 50 au is much less than that of the fiducial setup, while the total planetesimal mass of r_c = 200 au is approximately twice that of the fiducial setup.

For the central infall radius r_infall, the results for r_infall = 20 au (green) and r_infall = 80 au (orange) are presented. We note that the horizontal scale of the panel of r_infall in Figure 8 is different from the other panels in the same Figure. Due to constant σ, the widths of the gas surface density bumps (as well as the pressure bumps) increase with r_infall. The absolute values of negative η also increase with r_infall, because the initial gas surface density profile drops with increasing r. For r_infall = 20 au, planetesimal formation starts significantly earlier than the fiducial setup and produces a bit more planetesimal mass than the fiducial setup. For r_infall = 80 au, planetesimal formation starts much later than the fiducial setup, and produces much less planetesimal mass than the fiducial setup.

For the gas infall rate Ṁ_infall, the results for Ṁ_infall = 0.1 M_⊕ yr⁻¹ (green) and Ṁ_infall = 0.5 M_⊕ yr⁻¹ (orange) are presented. With the same initial disk conditions, the Gaussian-like gas surface density profile is higher and has a higher slope for larger Ṁ_infall than smaller Ṁ_infall, resulting in a broader region of negative η with higher absolute values. For Ṁ_infall = 0.1 M_⊕ yr⁻¹, no planetesimal is formed since ϵ_mid never reaches 1, though the condition for Q_p can be satisfied. For Ṁ_infall = 0.5 M_⊕ yr⁻¹, the results are quite similar to those of the simulation with M_disk = 0.02 M_⊙ except the total planetesimal mass. What the two systems have in common besides other parameters is the mass ratio between the infalling gas and the gaseous disk (Ṁ_infallt_infall)/M_disk. The difference in total planetesimal masses of the two systems is caused by different dust mass budgets.

For the infall duration t_infall, the results for t_infall = 20 kyr (green) and t_infall = 100 kyr (orange) are presented. Since the infall rate Ṁ_infall is constant, t_infall determines the total mass of infalling gas, leading to a higher gas surface density profile for larger t_infall. For t_infall = 20 kyr, the lengthy evolution time after the end of infall until t = 100 kyr flattens the gas surface density profile and makes the region with negative η negligible at this point. No planetesimal is formed during the evolution process since ϵ_mid never reaches 1, though the condition for Q_p can be satisfied. For t_infall = 100 kyr, the planetesimal formation process is much longer-lasting and more productive compared to the fiducial setup, though it starts slightly later.

For the relative infall width σ, the results for σ = 0.05 (green) and σ = 0.2 (orange) are presented. The gas surface density profiles present a similar trend as varying α: it is broader and lower for higher σ. For σ = 0.05, the negative part of η has higher absolute values than the fiducial setup, but spans a narrower range of r. The planetesimal formation process is longer-lasting than that of the fiducial setup, but less productive. For σ = 0.2, η is positive everywhere at this point, meaning that there is no pressure bump. Actually, a negligible pressure bump formed at 41 kyr and vanished at 74 kyr, but eventually no planetesimal is formed throughout the simulation since ϵ_mid never reaches 1. However, the condition for Q_p was still achieved in the simulation.

Fig. 9 Planetesimal mass evolution over 500 kyr for varying parameters compared with the fiducial setup. Each panel is labeled with the varied parameter, with green and orange curves representing decreased and increased parameter values in Table 1, respectively. Simulations that do not produce planetesimals are not displayed.

4.3 Sobol sampling analysis

We generated 4096 simulation setups, uniformly sampling the eight parameters within the value ranges shown in Table 1. A Sobol sequence (Sobol’ 1967) was used to achieve a quasirandom yet unclustered selection of parameter sets. Of the 4096 simulations, 653 successfully formed planetesimals. The histograms in Figure 10 illustrate the statistical distributions of these successful simulations across the eight parameters, and each parameter is divided into ten equal bins within its value range.

Among all parameters, the simulation outcomes are most sensitive to α, as evidenced by the steepest gradients in the first column of Figure 10 compared to the other parameters. At α ≈ 10⁻⁴, about 60% of the simulations result in planetesimal formation. This fraction drops rapidly as α increases. Simulations with α > 2.5 × 10⁻⁴ rarely produce planetesimals. The viscous diffusion timescale of the gas across the infall-induced pressure bump can be expressed as $$t_{\text{visc}}\simeq \frac{{\sigma }_r^2}{\upsilon }=\frac{{\sigma }^2{r}_{\text{infall}}^2}{\alpha c_{\text{s}}{H}_{\text{g}}}.$$(22)

A higher α leads to increased viscosity and a shorter gas diffusion timescale, which causes the pressure bump to smooth out more quickly. Conversely, with lower α, the pressure bump is retained longer, providing more time for dust to accumulate, grow, and trigger streaming instability.

The second most decisive parameter is v_frag. Simulations with v_frag < 200 cm s⁻¹ rarely form planetesimals. As v_frag increases, the maximum Stokes number of dust particles, determined by Eq. (20), also rises, leading to a higher average St for all particles. According to Eq. (9), particles with larger St are vertically distributed across a lower dust scale height H_d. Consequently, a higher average St increases the midplane dust density, ρ_d. This raises the midplane dust-to-gas ratio, ϵ_mid = ρ_d/ρ_g, facilitating the attainment of the critical condition for streaming instability, ϵ_mid = 1.

For M_disk and r_c, two parameters describing disk characteristics, the key to their effects on planetesimal formation is the initial gas mass within the infall region. According to Eq. (1), M_disk uniformly scales the magnitude of the initial gas surface density profile, while r_c determines how the initial gas mass is distributed across the disk radius r. A smaller M_disk reduces Σ_g,init, which leads to a steeper pressure bump under identical infall conditions. As r_c increases while M_disk remains fixed, the gas surface density decreases in the inner region of the disk and increases in the outer region, resulting in a flatter and more spatially extended profile. Since the range of r_infall tested is of the same order of magnitude as the minimum r_c, the infall region is located within the inner part of the disk. A larger r_c is conducive to forming the pressure bump, as it results in a lower initial gas surface density Σ_g,init at r_infall. If r_infall is located at hundreds of au, the effect of r_c will become more complex.

The infall location, r_infall, exhibits the weakest impact on the simulation outcomes among all parameters. According to Eq. (21), a lower orbital frequency Ω_K at a larger r_infall leads to a longer dust growth timescale, which hinders and postpones the onset of planetesimal formation. Nevertheless, the lower Σ_g,init at a larger r_infall promotes the formation of the pressure bump, similar to the effects of a smaller M_disk and a larger r_c. Consequently, the favorable and unfavorable factors for planetesimal formation at larger r_infall tend to offset each other. Although the histograms in Figure 10 show a slight preference for smaller r_infall, this parameter has a negligible impact on the likelihood of planetesimal formation and primarily influences its timing.

The other parameters describing infall properties, Ṁ_infall, t_infall, and σ, influence planetesimal formation by modulating the amount of infalling material and its distribution. The effects of Ṁ_infall and t_infall are primarily determined by their product, M_infall, which represents the total mass of the infalling gas. This is evidenced by the 2D histogram of Ṁ_infall and t_infall in Figure 10. Since only two of these three quantities are independent, we selected the two that are mostly measured in recent observational studies of infall streamers (Alves et al. 2020; Pineda et al. 2020; Cacciapuoti et al. 2024). A larger Ṁ_infall or t_infall leads to a greater M_infall, which not only enhances the prominence of the pressure bump, but also provides more dust material to build planetesimals. The relative width of the infall region, σ, determines how concentrated M_infall is distributed over r. A lower σ creates a steeper but narrower pressure bump. While the steepness aids in triggering streaming instability, the narrower width limits the amount of material available for planetesimal formation, which explains the intersection of the two curves in the panel of σ in Figure 9.

In summary, planetesimal formation in an infall-induced pressure bump, with r_infall located within a few tens of au, is more likely to occur with lower values of α, M_disk, and σ, and higher values of v_frag, r_c, Ṁ_infall, and t_infall. The most decisive two factors are α and v_frag, as planetesimal formation becomes highly improbable when α > 2.5 × 10⁻⁴ and v_frag < 200 cm s⁻¹. The ideal scenario is a low-viscosity, low-mass, but spatially extended disk accreting a massive yet narrow infall streamer.

Fig. 10 Histograms of 653 simulations where planetesimal formation occurs, out of 4096 Sobol-sampled simulations. “Counts” refers to the number of simulations. The 1D histograms along the diagonal show the distribution over each of the eight parameters individually. The 2D heatmap histograms display the joint distributions over pairs of parameters, with gray cells indicating a count of 0.

5 Discussions

5.1 Timescale ratios

In order to reach the critical conditions for planetesimal formation, a protoplanetary disk undergoing infall must satisfy the following conditions: (1) forming a pressure bump prominent enough to trap and accumulate dust; and (2) ensuring that the pressure bump persists long enough to allow dust to grow until the onset of streaming instability. Each of these two conditions can be characterized by a timescale ratio, which is calculated by the initial parameter setups of the simulations.

The presence of a pressure bump is signified by a super- Keplerian region with positive pressure gradient and negative η values, located on the inner side of the bump. According to Eq. (11), the minimum value of η (i.e., the largest absolute value for negative η) corresponds to the fastest outward drift speed of dust, reflecting the pressure bump’s ability to accumulate and trap dust. The dust drift timescale at this site is formulated as $$t_{\text{drift}}=\frac{r}{{\upsilon }_{\text{d},\text{max}}}=-\frac{1+{\text{St}}^2}{2{\eta }_{\min }\text{St}{\mathrm{\Omega }}_{\text{K}}}.$$(23)

Given the disk parameters M_disk and r_c, along with the infall parameters Ṁ_infall, t_infall, r_infall, and σ, the gas surface density profile after the infall process, assuming no viscosity, is calculated as $${\text{Σ}}_{\text{g}}={\text{Σ}}_{\text{g},\text{\hspace{0.17em}init\hspace{0.17em}}}+S_{\text{g},\text{\hspace{0.17em}infall\hspace{0.17em}}}{t}_{\text{infall}}.$$(24)

Using the relations from Eqs. (3) to (5), the analytical radial profile of η can be derived, allowing for the calculation of its minimum value η_min. The Stokes number at the fragmentation limit at this location, St_frag, given by Eq. (20), is used in the calculation. If no fragmentation limit is present, St_frag = 1 is assumed, representing the largest particles participating in streaming instability, as described in Section 2.3. Given the dust growth timescale in Eq. (21), the timescale ratio of drift to growth is expressed as $$\frac{t_{\text{drift\hspace{0.17em}}}}{t_{\text{grow\hspace{0.17em}}}}=-\frac{\left(1+{\text{St}}_{\text{frag}}^2\right)ϵ}{2{\eta }_{\min }{\text{St}}_{\text{frag\hspace{0.17em}}}},$$(25)

where a constant initial global dust-to-gas ratio ϵ = ϵ_init = 0.01 is assumed, since ϵ evolves over time and varies with r. This quantity is influenced by all the eight major parameters.

The lifetime of a pressure bump is measured by the gas viscous diffusion timescale across the bump, given by Eq. (22). The timescale ratio of gas viscous diffusion to dust growth is formulated as $$\frac{t_{\text{visc}}}{t_{\text{grow\hspace{0.17em}}}}=\frac{{\sigma }^2{r}_{\text{infall\hspace{0.17em}}}^2{\mathrm{\Omega }}_{\text{K}}^2ϵ}{\alpha c_{\text{s}}^2}\propto \frac{{\sigma }^2}{\alpha r_{\text{infall\hspace{0.17em}}}^{1/2}},$$(26)

where it is considered that c_s and Ω_K vary only with r_infall. Comparing t_visc to t_grow reveals the extent to which dust particles can grow before the pressure bump dissipates. In contrast, comparing t_drift to t_grow serves more for mathematical clarity than physical necessity, as t_drift alone effectively measures the pressure bump’s ability to trap dust.

Figure 11 displays the distribution of the two timescale ratios in logarithmic scale, calculated for our simulations. The blue and orange markers represent the Sobol-sampled simulations, excluding cases where the pressure bump fails to form (where t_drift is negative). The two timescale ratios are not independent: a lower α increases St_frag, leading to a lower t_drift /t_grow and a higher t_visc /t_grow ; a lower σ steepens the pressure gradient (increasing |η_min|), which reduces both t_drift/t_grow and t_visc/t_grow. A distinct boundary separates simulations that form planetesimals (blue) from those that do not (orange). For a given t_visc /t_grow, simulations with faster outward dust drift in the pressure bump (lower t_drift/t_grow) are more likely to produce planetesimals. Likewise, for a given t_drift/t_grow, planetesimal formation is favored by longer diffusion timescales, indicating a longer-lived pressure bump relative to dust growth timescales. A logistic regression fit gives the boundary as: $${\log }_{10}\left(\frac{t_{\text{visc\hspace{0.17em}}}}{t_{\text{grow\hspace{0.17em}}}}\right)>1.00{\log }_{10}\left(\frac{t_{\text{drift\hspace{0.17em}}}}{t_{\text{grow\hspace{0.17em}}}}\right)+1.73,$$(27) $$t_{\text{visc\hspace{0.17em}}}\gtrsim 53.55t_{\text{drift\hspace{0.17em}}},$$(28)

which defines the conditions for planetesimal formation, highlighting the need for a long-lasting pressure bump that is sufficiently prominent to facilitate rapid dust accumulation.

The distribution of blue and orange markers leaves an empty region in the lower-left corner of Figure 11. To investigate whether there is a lower limit for t_visc /t_grow within the planetesimal-forming domain, we selected simulations in the range −0.75 < log₁₀(t_drift/t_grow) < −0.5 and 1.05 < log₁₀ (t_visc /t_grow) < 1.45, where most simulations form planetesimals. We then reduced their initial global dust-to-gas ratio to ϵ_init = 10^−2.25, filling the previously blank area. These additional simulations are shown as green and red markers in Figure 11. Although many of the red crosses lie to the upper left of the boundary line given by Eq. (27), they do not result in planetesimal formation. This suggests that if the pressure bump’s lifetime is insufficient relative to the dust growth timescale, planetesimals cannot form in time, even if dust accumulation via outward drift in the pressure bump is sufficiently rapid. We approximate this boundary as: $$t_{\text{visc\hspace{0.17em}}}\gtrsim 10t_{\text{grow\hspace{0.17em}}}.$$(29)

Combining this with Eq. (28) provides the critical conditions for planetesimal formation in an infall-induced pressure bump.

Fig. 11 Distribution of timescale ratios tdrift/tgrow and tvisc /tgrow. The plus symbols (“+”) represent simulations in which planetesimals form, while the cross symbols (“x”) indicate simulations with no planetesimal formation. Blue and orange markers show the Sobol-sampled simulations, while red and green markers represent a subset of 30 selected simulations from these Sobol-sampled simulations, with a reduced initial global dust-to-gas ratio (ϵinit = 10−2.25). The dashed gray line denotes the fitted boundary of the planetesimal-forming domain.

5.2 Comparison of contributions from primordial dust and infalling dust

In our models, the infalling streamer shares the same global dust-to-gas ratio as the initial protoplanetary disk, meaning that the solid material forming planetesimals originates from both the dust in the original disk and the infalling dust. To determine whether planetesimals can form with only one dust source, and which source contributes more, we simulated two extreme scenarios: one with a dust-free initial disk and another with a dust-free infall streamer. This experiment was conducted with α = {1, 1.5, 2} × 10⁻⁴, while keeping all other parameters at their fiducial values.

The results are summarized in Table 3. For α = 10⁻⁴, planetesimal formation occurs in both scenarios. The combined planetesimal mass in these two extreme scenarios is less than that of the simulation with both dust sources present. For α = 1.5 × 10⁻⁴, planetesimal formation occurs only in the simulation with an initially dust-free disk, where the planetesimal mass formed from infalling dust alone is less than in the simulation with both dust sources. For α = 2 × 10⁻⁴, planetesimal formation does not occur in either extreme scenario. Therefore, planetesimal formation with only one dust source is possible when conditions such as a low α value are particularly favorable, with the infalling dust contributing more significantly than the primordial dust in the disk.

Once an infall-induced pressure bump forms, all incoming infalling dust is trapped by the bump, whereas only a portion of the dust initially in the disk can be captured. Since the mass of the infalling dust is comparable to the total initial dust mass in the disk, the infalling dust becomes the primary contributor to planetesimal formation. However, the total amount of planetesimals formed depends on multiple factors beyond just the amount of available dust. The standard scenario, which includes both dust sources, not only traps more dust in the pressure bump but also reaches the critical conditions for streaming instability earlier than the two extreme scenarios with only one dust source. This explains why the standard scenario produces more planetesimal mass than the combined total of the two extreme scenarios.

In all the simulations above, infall begins at t = 0. We also tested scenarios where infall starts later. Since dust growth begins from sub-micron-sized particles, a delayed infall means that dust grains initially in the disk have already grown to larger sizes when the pressure bump forms, resulting in faster dust drift and trapping in the pressure bump. However, if infall starts too late, a significant portion of the primordial dust drifts inward to the inner disk in the absence of a dust trap, reducing the amount of dust available to be trapped by the pressure bump. We tested starting infall at 50 kyr and 100 kyr on the fiducial setup. In the former case, the total planetesimal mass increases to 11.57 M_⊕ due to the faster accumulation of the primordial dust. In the latter, it drops to 9.71 M_⊕ due to the reduced amount of trapped primordial dust. This effect is minor since the infalling dust plays a dominant role, and the timing of infall only affects the distribution of the dust initially present in the disk.

5.3 Effects of other parameters

Apart from the eight major parameters we previously discussed, two other parameters also play roles in shaping the planetesimal formation process: the small-scale diffusion parameter δ (Eq. (15)) and the planetesimal formation efficiency ζ (Eq. (17)). These parameters were excluded from our parameter space because they are not directly relevant to the large-scale evolution of gas and dust in the disk.

The small-scale diffusion parameter δ is defined as the ratio of radial particle diffusivity to c_sH_g (Schreiber & Klahr 2018), analogous to the definition of α in terms of gas viscosity. According to Eq. (15), a decrease in δ leads to a lower Toomre-like parameter Q_p, which facilitates reaching the critical condition for streaming instability, Q_p < 1. For all simulations presented above, we used δ = 10⁻⁴, the most conservative value for planetesimal formation according to the measurements by Schreiber & Klahr (2018). In contrast, recent work by Lau et al. (2024), which employs the same criterion for planetesimal formation, adopts δ = 10⁻⁵. However, our simulation results indicate that the other criterion, ϵ_mid ≥ 1, is always more restrictive than Q_p < 1. Thus, reducing δ does not affect the triggering of streaming instability; its only effect in this context is to increase the activation function 𝒫_pf through the reduced Q_p, as shown in Eq. (18), ultimately influencing the total planetesimal mass formed, as described by Eq. (17). We tested δ = 10⁻⁵ on the fiducial setup and found that the start and end times of planetesimal formation were identical to those of the fiducial setup, with a slight increase in total planetesimal mass, from 9.97 M_⊕ to 9.99 M_⊕. Besides determining Q_p, δ also plays a key role in defining the initial mass function of planetesimals, as given by Gerbig & Li (2023). This aspect, however, is not addressed in this work, as we do not convert planetesimal surface density into planetary bodies.

The planetesimal formation efficiency ζ also has no influence on the triggering of streaming instability. As shown in Eq. (17), ζ determines the mass of planetesimals formed per unit time, with higher values of ζ resulting in faster planetesimal formation. Drążkowska et al. (2016), who propose this prescription for planetesimal transformation, demonstrate that a range of ζ values between 10⁻⁶ and 10⁻² is reasonable. In all simulations presented so far, we adopted ζ = 10⁻³. We tested ζ = 10⁻² on the fiducial setup. The total planetesimal mass formed is 24.67 M_⊕, significantly higher than in the fiducial simulation. While the onset time of planetesimal formation remains the same as in the fiducial setup, the end time occurs earlier, shortening the planetesimal formation duration from 162 kyr to 119 kyr. Since planetesimal formation ceases when ϵ_mid drops below one, a faster rate of dust consumption shortens the formation duration. Within a given time frame, a quicker conversion rate results in a greater total mass of planetesimals. The overall effect increases the total planetesimal mass, though not as significantly as the amplification of ζ.

The radial resolution is a numerical factor that can influence simulation results. To assess its impact, we increased the number of refined cells in the infall region from 60 to 80 and 100 in the fiducial simulation. The total mass of planetesimals formed increased slightly from 9.97 M_⊕ to 10.80 M_⊕ and 11.32 M_⊕, respectively. Although the results are not fully convergent with respect to radial resolution, this minor variation does not impact our conclusions.

5.4 Caveats

5.4.1 Rossby wave instability

The Rossby wave instability (RWI) can be triggered in a narrow, axisymmetric pressure bump by strong pressure gradients and non-Keplerian radial shear (Lovelace et al. 1999; Li et al. 2000). RWI generates anticyclonic vortices, which can disrupt the axisymmetric ring-like pressure bump. However, these vortices create azimuthal pressure bumps that efficiently trap dust (Meheut et al. 2012), where the dust-to-gas ratio may exceed that of an axisymmetric pressure bump, potentially triggering streaming instability.

To assess the stability of the infall-induced pressure bumps in this study against RWI, we applied the criterion for RWI proposed by Chang et al. (2023), $${\kappa }^2+N_r^2<0.6{\mathrm{\Omega }}_{\text{K}}^2,$$(30)

where κ is the epicyclic frequency, and N_r is the radial buoyancy frequency.

In our fiducial simulation, an RWI-unstable region forms between 40.5 and 42.9 au, persisting from 14 to 51 kyr. Assuming no diffusion (calculating Σ_g using Eq. (24)), 1788 of the 4096 Sobol-sampled simulations remain stable against RWI, and 215 of the 653 simulations where planetesimal formation occurs are stable. This is a conservative estimate, as diffusion reduces the magnitude of the pressure gradient, lowering the likelihood of RWI.

Simulations have demonstrated that RWI vortices driven by infall are possible (Bae et al. 2015; Kuznetsova et al. 2022). However, the impact of infall-induced RWI on planetesimal formation – whether positive or negative – remains uncertain and requires further investigation through 2D or 3D simulations.

5.4.2 Planetary evolution

In our simulations, the planetesimal surface density profile and total planetesimal mass cease evolving once the conditions for streaming instability are no longer satisfied. However, a planetesimal belt with a width of only a few au and a total mass of just a few Earth masses is unlikely the ultimate fate of the system, as planetary evolution and its feedback on the protoplanetary disk have not been considered.

Since the planetesimals are formed in a dust-rich environment, they are expected to grow efficiently through pebble accretion (Ormel & Klahr 2010; Lambrechts & Johansen 2012). Pebble accretion might lead to a continuous increase in the total planetesimal mass even after the formation phase ends, accelerate the depletion of dust mass, and cause an earlier cessation of planetesimal formation due to rapid dust consumption.

Gas accretion begins when a planetary embryo reaches the pebble isolation mass (Lambrechts et al. 2014), transforming a fraction of the gas in the disk into planetary mass. A sufficiently massive planet can carve a gap in the disk (Kanagawa et al. 2017; Duffell 2020), creating a new pressure bump outside the gap. This might lead to the formation of a new generation of planetesimals, resulting in sequential planetesimal formation, as explored by Lau et al. (2024).

The aerodynamic drag (Adachi et al. 1976) and planet-disk interaction torques (Goldreich & Tremaine 1979; Tanaka et al. 2002) generally cause planets to lose angular momentum and migrate inward. The pressure bump induced by infall might serve as a migration trap eliminating such torques and retaining planets (Coleman & Nelson 2016; Morbidelli 2020). The growth of planets and their interactions with the disk following planetesimal formation in an infall-induced pressure bump merit further investigation through N-body simulations.

6 Conclusions

This work demonstrates the feasibility of planetesimal formation from dust in a pressure bump created by the late-stage infall of interstellar medium onto a protoplanetary disk. We used DustPyto simulate gas evolution, dust transport and growth, and the transformation of dust into planetesimals triggered by streaming instability.

The Gaussian-distributed infalling gas generates a local maximum in the radial profile of gas surface density, leading to a local maximum in midplane pressure slightly inward of the gas surface density peak. This pressure bump halts the inward drifting dust from the outer disk and traps the infalling dust, enhancing the local dust-to-gas ratio and accelerating dust growth. When the critical conditions for streaming instability, governed by the midplane dust-to-gas ratio, are satisfied, planetesimals begin to form continuously, converting the dust surface density into planetesimal surface density. As the pressure bump is gradually smoothed out by viscous diffusion of gas, planetesimal formation ceases when the dust trap vanishes and the peak midplane dust-to-gas ratio drops below the threshold for streaming instability. Ultimately, the viscous spreading of the bump results in the outside-in formation of a planetesimal belt inside the central infall location.

A large number of simulations were conducted to investigate the effects of various parameters related to disk and infall properties. These parameters influence whether planetesimal formation occurs, the timing of formation, and the total mass of planetesimals produced. Planetesimal formation is favored by conditions where a massive, narrowly distributed infall of interstellar medium occurs onto a disk with low viscosity, low mass, extended surface density distribution, and high dust fragmentation velocity. Formation is highly unlikely when α > 2.5 × 10⁻⁴ due to rapid gas diffusion, or when v_frag < 200 cm s⁻¹ due to the fragility of dust particles. The critical conditions for planetesimal formation within an infall-induced pressure bump can be characterized by specific timescale relationships: t_visc ≳ 53.55 t_drift, t_visc ≳ 10 t_grow. These timescales are directly determined by the system’s parameter setup. These quantitative conditions encapsulate the essential qualitative requirements for planetesimal formation: the creation of a prominent pressure bump that efficiently traps dust and its prolonged persistence. The primary source of material for planetesimals is the infalling dust; however, under favorable conditions, planetesimal formation can still occur even if the infalling material is solely gas.

In this study, we focused on the conditions and processes leading to planetesimal formation, but the subsequent stages of planetary evolution remain unexplored. Future work will extend our model to include the growth of planets through pebble and gas accretion, the dynamical interactions between planets, and their interactions with the gaseous disk. To achieve this, we need to couple our current dust and gas evolution code with an N-body simulation framework, allowing us to comprehensively model the evolution of a planetary system starting from the initial formation of planetesimals within an infall-induced pressure bump.

Data availability

Movie associated to Section 4.1 is available at https://www.aanda.org

Acknowledgements

We acknowledge funding from the European Union under the European Union’s Horizon Europe Research and Innovation Programme 101124282 (EARLYBIRD) and 101040037 (PLANETOIDS) and funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under grant 325594231, and Germany’s Excellence Strategy – EXC-2094 – 390783311. Views and opinions expressed are, however, those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them.

References

All Tables

All Figures

	Fig. 1 Illustration of the infall model, showing profiles of gas surface density (blue curves) and pressure (orange curves). The dashed curves represent the initial profiles, while the solid curves depict the profiles after infall. The horizontal blue and orange arrows indicate the directions of gas and dust radial velocities, respectively.
In the text

	Fig. 2 Smooth activation functions in terms of Qp with smoothness parameter n = 0.2, centered at Qp= 0.75. Gravitational collapse occurs when Qp < 1.
In the text

	Fig. 3 Dust distribution of the fiducial simulation at the start of planetesimal formation. The dust surface density σd is independent of the mass grid. The green contour indicates the drift limit, while the blue contour marks the fragmentation limit. White contours represent particle masses corresponding to Stokes numbers St = 1 (solid) and St = {10−1, 10−2, 10−3 } (dotted), respectively.
In the text

Fig. 4 Midplane dust-to-gas ratio ϵmid, Toomre-like value Qp, and pressure gradient parameter η around the infall region of the fiducial simulation at the start of planetesimal formation. The vertical dotted gray line indicates rinfall = 40 au. The horizontal dashed gray line represents the critical conditions for streaming instability, where ϵmid = 1 and Qp = 1, as well as the condition for the presence of a pressure bump, η = 0.

In the text

	Fig. 5 Radial surface density distributions of gas, dust, and planetesimals, as well as the midplane dust-to-gas ratio ϵmid for the fiducial simulation at the start (top) and end (bottom) of planetesimal formation. The vertical dotted gray line indicates rinfall = 40 au. The horizontal dashed gray lines represent ϵmid = 1.
In the text

	Fig. 6 Mass evolution of the gas, dust, and planetesimals over 500 kyr in the fiducial simulation. Dotted lines indicate the start and end of planetesimal formation.
In the text

	Fig. 7 Gas surface density distributions at t = 100 kyr for varying parameters compared with the fiducial setup. Each panel is labeled with the varied parameter, with green and orange curves representing decreased and increased parameter values in Table 1, respectively.
In the text

	Fig. 8 Pressure gradient parameter η at t = 100 kyr for varying parameters compared with the fiducial setup. Each panel is labeled with the varied parameter, with green and orange curves representing decreased and increased parameter values in Table 1, respectively. The horizontal dashed gray lines represent η = 0.
In the text

	Fig. 9 Planetesimal mass evolution over 500 kyr for varying parameters compared with the fiducial setup. Each panel is labeled with the varied parameter, with green and orange curves representing decreased and increased parameter values in Table 1, respectively. Simulations that do not produce planetesimals are not displayed.
In the text

	Fig. 10 Histograms of 653 simulations where planetesimal formation occurs, out of 4096 Sobol-sampled simulations. “Counts” refers to the number of simulations. The 1D histograms along the diagonal show the distribution over each of the eight parameters individually. The 2D heatmap histograms display the joint distributions over pairs of parameters, with gray cells indicating a count of 0.
In the text

Fig. 11 Distribution of timescale ratios tdrift/tgrow and tvisc /tgrow. The plus symbols (“+”) represent simulations in which planetesimals form, while the cross symbols (“x”) indicate simulations with no planetesimal formation. Blue and orange markers show the Sobol-sampled simulations, while red and green markers represent a subset of 30 selected simulations from these Sobol-sampled simulations, with a reduced initial global dust-to-gas ratio (ϵinit = 10−2.25). The dashed gray line denotes the fitted boundary of the planetesimal-forming domain.

In the text