© The Authors 2024

1 Introduction

Observations of protoplanetary disks have revealed substructures in dust profiles at distances greater than 10 au (e.g., ALMA Partnership 2015; Andrews et al. 2018). An unbiased proto- planetary disk survey in the Taurus star-forming region, where approximately 75% of solar-mass stars have disks (Luhman et al. 2009), exhibits a substructure occurrence rate as high as 40% (Long et al. 2018). The most common type of dust substructures are annular depletions and enhancements in the continuum emissions, which are referred to as dust rings and gaps.

Several mechanisms have been proposed to explain the dust rings and gaps, such as various types of instabilities, processing of dust at the snow lines, magnetohydrodynamic effects, and pressure maxima in a radial pressure profile (Bae et al. 2023, and references therein). In addition to the preceding mechanisms, a widely accepted formation channel is by planets carving gas gaps with masses typically ≳15 M_⊕ (Earth masses). Hereafter we refer to this mechanism as the gas-gap mechanism (e.g., Paardekooper & Mellema 2006). If an observed dust gap at ≳10 au is caused by an unseen planet, the planet mass can be estimated from the results of the disk–planet interaction simulations (Zhang et al. 2018; Lodato et al. 2019; Wang et al. 2021). The inferred masses of putative planets are distributed in a range of a few Earth masses to ∼10 Jupiter masses, ∼70% of which have >0.1 Jupiter mass (Bae et al. 2023). Considering the fraction of disks with substructures, these estimates suggest that the occurrence fraction of planets with masses exceeding 0.1 Jupiter mass in wide orbits is approximately 20%. This value appears to be in tension with the low occurrence rate of cold gas giants as suggested by the current observed period-mass diagram of exoplanets (<10%; Fernandes et al. 2019; Fulton et al. 2021) and predicted occurrence rates in population synthesis models (Mordasini et al. 2018; Emsenhuber et al. 2021). The current period and mass distribution of exoplanets could be reproduced if these putative planets undergo large-scale inward migration (Lodato et al. 2019; Mulders et al. 2021; van der Marel & Mulders 2021). However, the feasibility of this scenario could be low due to inefficient type II migration in low-viscosity disks (Ndugu et al. 2019; Müller-Horn et al. 2022; Tzouvanou et al. 2023).

Dust substructures can be created by low-mass, no-gas-gapopening planets (≲10M_⊕) in disks. A dust gap forms due to the gravitational interaction between the planet and the dust (Muto & Inutsuka 2009; Dipierro et al. 2016; Dipierro & Laibe 2017). In our previous work, Kuwahara et al. (2022) (hereafter Paper I), we showed that the gas flows driven by low-mass planets can create dust substructures in disks with low turbulent viscosity at the disk midplane (hereafter referred to as the gas-flow mechanism). A low-mass protoplanet (typically ~0.1–10 M_⊕) embedded in a disk induces a three-dimensional (3D) gas flow (e.g., Ormel et al. 2015; Fung et al. 2015; Kuwahara et al. 2019). If the disk midplane is weakly turbulent, as suggested by recent studies (Villenave et al. 2022; Jiang et al. 2024), the radially-outward outflow of the gas generates a congestion of dust outside the planetary orbit, because the radially-inward outflow blows dust away from the planetary orbit. This leads to the formation of a dust ring outside the planetary orbit, and a gap interior to it. This mechanism thus differs from the dust substructures generated by carving gas gaps, as occurs around higher mass planets. The dust ring and gap formation by low-mass, no-gas-gap-opening planets could therefore reconcile the frequently observed dust gaps seen in disks that have no corresponding gas gaps (Zhang et al. 2021; Jiang et al. 2022). Moreover, it may be a fresh perspective on the proposed large fraction of low-mass planets (≲10 M_⊕) at wide orbits (≳ 10³ days) inferred from a population synthesis model (Drazkowska et al. 2023).

Paper I assumed a steady state for simplicity and computed the dust surface density perturbed by the planet-induced gas flow. However, the validity of the steady-state assumption is nontrivial. The 3D structure of the gas flow has a complex dependence on different parameters such as the planetary mass and the pressure gradient of the disk gas (Ormel et al. 2015; Kurokawa & Tanigawa 2018), which in turn regulates how dust piles up outside the orbit of the planet and is depleted interior to it. It is therefore still unclear how the profiles of dust rings and gaps vary with time and compare with observed disks whose ages are typically a few Myr (Haisch et al. 2001).

In this study, we investigate the time evolution of the dust surface density perturbed by the planet-induced gas flow (Sect. 3). We extend the parameter space and conduct a more comprehensive investigation than in Paper I. In addition, we introduce semi-analytic models describing the properties of the dust rings and gaps, such as the widths and the depths based on the results of numerical simulations (Sect. 4.4). These approaches allow us to efficiently explore the disk parameter space where low- mass planets can create dust gaps and rings that are comparable in magnitude to those observed in young disks. In Sect. 5 we discuss the implications for planet formation and observations of protoplanetary disks, showing that up to ∼65% (∼15%) of the observed dust gaps (rings) could be caused by the gas flow induced by low-mass planets at wide-orbits. We conclude in Sect. 6. Our key results can be seen in Eqs. (26), (28), and (33), which are the semi-analytic models describing the widths of the dust ring and gap and the depth of the dust gap. We compare these models with the observational data in Figs. 19 and 20 in Sect. 5.

2 Numerical methods

In Paper I we investigated the dust substructure formation by the gas-flow mechanism in three steps: (1) we first performed 3D hydrodynamical simulations of the gas flow around an embedded planet and obtained the gas flow velocity field. (2) By the post-processing these simulations, we calculated the radial drift velocity of dust perturbed by the gas flow, where the obtained gas velocity field was used to compute the dust motion. (3) Finally, assuming a steady state, we computed the dust surface density by incorporating the obtained perturbed radial drift velocity of dust into a one-dimensional (1D) advection-diffusion equation (see also Fig. 1 of Paper I).

In this study, we followed the same procedures described above, but we investigated the time-dependent dust surface density in step 3 described above. The following sections summarize our numerical approach.

2.1 Non-dimensionalization

As in Paper I, in our simulations the length, times, velocities, and densities are normalized by the disk gas scale height, H, the reciprocal of the orbital frequency, Ω⁻¹, isothermal sound speed, c_s, and the unperturbed gas density at the location of the planet, ρ_∞.

In this dimensionless unit system, we defined the dimensionless thermal mass of the planet, $$m\equiv \frac{R_{\text{Bondi\hspace{0.17em}}}}{H}=\frac{M_{\text{p}}}{M_{\text{th}}},$$(1)

where $R_{\text{Bondi\hspace{0.17em}}}=G{M}_{\text{p}}/c_{\text{s}}^2$ is the Bondi radius, G is the gravitational constant, M_p is the mass of the planet, M_th = M_*h³ is the thermal mass, M_* is the stellar mass, and h is the aspect ratio of the disk. The Hill radius is given by R_Hill = (m/3)^1/3 H. In Paper I, we considered three planetary masses: m = 0.03, 0.1, and 0.3. In this study, we considered eight planetary masses ranging from m = 0.03 to 0.5 (Table 1), which corresponds to planets with M_p ≃ 0.2–3.3 M_⊕ orbiting a solar-mass star at 1 au (M_p ≃ 3.9–66 M_⊕ at 50 au; Eq. (A.12)). Throughout the paper, when we convert the dimensionless quantities into dimensional ones, we considered the typical steady accretion disk model with a constant turbulence strength (Shakura & Sunyaev 1973), including viscous heating due to the accretion of the gas and irradiation heating from the central star (Appendix A; Ida et al. 2016).

Our planet revolves with the Keplerian speed, v_K, on a fixed circular orbit. Because the disk gas rotates with the sub- Keplerian speed due to the global pressure gradient, planet experiences the headwind of the gas. We defined the Mach number of the headwind as $${ℳ}_{\text{hw}}\equiv -\frac{h}{2}\left(\frac{\text{\hspace{0.17em}d\hspace{0.17em}ln\hspace{0.17em}}p}{\text{\hspace{0.17em}d\hspace{0.17em}ln\hspace{0.17em}}r}\right),$$(2)

where p is the pressure. We considered three Mach numbers: ℳ_hw = 0.01, 0.03, and 0.1 (Fig. A.1).

The global pressure gradient of the disk gas causes the radial drift of dust. The unperturbed drift velocity is given by (Weidenschilling 1977; Nakagawa et al. 1986) $${\upsilon }_{\text{drift}}=-\frac{2\text{St}}{1+{\text{St}}^2}{ℳ}_{\text{hw}},$$(3)

where $$\text{St}=t_{\text{stop\hspace{0.17em}}}\text{Ω}$$(4)

is the Stokes number of dust and t_stop is the stopping time of dust. Because Paper I found that the apparent dust ring and gap form when St ≲ 10⁻², we considered St = 10⁻⁴–10⁻², which corresponds to ∼0.37–37 mm sized dust grains at 1 au (∼ 7.6 × 10⁻³–0.76 mm at 50 au; Appendix A).

2.2 Hydrodynamical simulations

Assuming a compressible, inviscid, non-self-gravitating sub- Keplerian gas disk with the vertical stratification due to the stellar gravity, we performed 3D nonisothermal hydrodynamical simulations using Athena++ code¹ (Stone et al. 2020). Our methods of hydrodynamical simulations are the same as described in Paper I, but this study handles a broader and more detailed parameter space compared to Paper I in terms of the planetary mass (Table 1). Our hydrodynamical simulations were performed in the local frame co-rotating with the planet (see also Fig. 1 of Paper I). Radiative cooling was implemented by using the so-called β cooling model, where the radiative cooling occurs on a finite timescale, β (Gammie 2001). Following Kurokawa & Tanigawa (2018), we set the dimensionless cooling time as β = (m/0.1)². We simulated the gas flow for at least 10² Keplerian orbits, where the flow field seems to have reached a steady state (see Sect. 2.3 of Paper I, for details).

2.3 Calculations of the radial radial drift velocity of dust perturbed by the gas flow

The radial drift velocity of dust perturbed by the planet-induced gas flow was calculated by the same method as Paper I. (1) We first numerically integrated the equation of motion of dust in a local domain co-rotating with the planet (the local Cartesian coordinates (x, y, z) centered at the planet), in which we used the gas velocity obtained from the hydrodynamical simulation to calculate the gas drag force acting on dust (Kuwahara & Kurokawa 2020). Hereafter we denote the x-, y-, and z-directions as the radially outward, azimuthal and vertical directions to the disk, respectively. (2) We then sampled the positions and the velocities of dust at fixed small time intervals in the local domain of orbital integration of dust, obtaining in this way the spatial distribution of dust. (3) We assumed the uniform and Gaussian distributions of dust in the azimuthal and vertical directions outside the local domain of orbital integration of dust, in which dust has the unperturbed steady-state drift velocity, v_drifte_x . (4) Finally, we computed the radial drift velocity of dust perturbed by the planet-induced gas flow, ⟨v_d⟩, by averaging the x-component of the dust velocity in the vertical and full azimuthal directions in a disk (see Sects. 2.4–2.5 of Paper I, for details).

2.4 Dust surface density calculation

We computed the dust surface density by incorporating the perturbed radial drift velocity of dust into a 1D advection-diffusion equation: $$\frac{\partial {\text{Σ}}_{\text{d}}}{\partial t}+\frac{\partial }{\partial x}\left(\langle v_{\text{d}}\rangle {\text{Σ}}_{\text{d}}-\mathcal{D}\frac{\partial {\text{Σ}}_{\text{d}}}{\partial x}\right)=0.$$(5)

Here Σ_d is the dust surface density, 𝒟 = α_diff/(1 + St²) is the diffusion coefficient for the dust (Youdin & Lithwick 2007), and α_diff is a dimensionless turbulent parameter describing turbulent diffusion of dust. Because Paper I found that the dust rings and gaps do not appear when α_diff ≳ 10⁻³, in this study we only assumed α_diff = 10⁻⁴ (hereafter referred to as the moderate-turbulence case) and 10⁻⁵ (hereafter referred to as the low-turbulence case). In Eq. (5), we neglect the effect of the disk curvature by focusing on a radial range sufficiently narrow compared to the orbital radius of the planet. We did not consider the backreaction of dust on gas.

While Paper I assumed a steady state in Eq. (5), we computed the time evolution of the dust surface density in this study. We assumed that the gas surface density is constant for simplicity, so that Eq. (5) does not contain the gas surface density. We integrated Eq. (5) using a finite-volume method. The size of the calculation domain of dust surface density simulation was x є [x_in, x_out]. We set x_in = −100 and x_out = 100. We used a fixed spatial interval Δx = 0.01. A zero-diffusive flux condition was adopted at x = x_in. At x = x_out, we set the constant advective flux, v_driftΣ_d,0, where Σ_d,0 = 1. The time step was calculated as $$\text{Δ}t=\text{CFL}\times \min \left(\frac{1}{\underset{{}_i}{\max }\left(\left|{\langle v_{\text{d}}\rangle }_i\right|/\text{Δ}x\right)},\frac{{(\text{Δ}x)}^2}{\mathcal{D}}\right),$$(6)

where the Courant–Friedrichs–Lewy (CFL) number was set to CFL = 0.5 and 〈v_d〉_i is the dust velocity at the i-th grid.

2.5 Analytic formulae for analyzing the numerical results

In the following sections (Sects. 2.5.1–2.5.4), we introduce the analytic formulae which will be used to analyze the results of numerical simulations.

2.5.1 Drift and diffusion timescales of dust

We identified two key timescales: the drift timescale of dust $$t_{\text{drift\hspace{0.17em}}}=\frac{ℒ}{\left|{\upsilon }_{\text{drift\hspace{0.17em}}}\right|}\simeq \frac{ℒ}{2\text{St}{ℳ}_{\text{hw}}}\simeq 1.67\times {10}^4{\left(\frac{\text{St}}{{10}^{-3}}\right)}^{-1}{\left(\frac{{ℳ}_{\text{hw}}}{0.03}\right)}^{-1}\left(\frac{ℒ}{H}\right)$$(7)

and the diffusion timescale $$t_{\text{diff\hspace{0.17em}}}=\frac{{ℒ}^2}{\mathcal{D}}\simeq \frac{{ℒ}^2}{{\alpha }_{\text{diff\hspace{0.17em}}}}={10}^4{\left(\frac{{\alpha }_{\text{diff\hspace{0.17em}}}}{{10}^{-4}}\right)}^{-1}{\left(\frac{ℒ}{H}\right)}^2.$$(8)

Here ℒ is the characteristic length and we assumed 1 + St² ≃1 (St ≪ 1). The drift timescale coincides with the diffusion timescale when $$ℒ=\frac{{\alpha }_{\text{diff}}}{2\text{St}{ℳ}_{\text{hw}}}\equiv {ℒ}_{\text{eq}}\simeq 1.67\left(\frac{{\alpha }_{\text{diff}}}{{10}^{-4}}\right){\left(\frac{\text{St}}{{10}^{-3}}\right)}^{-1}{\left(\frac{{ℳ}_{\text{hw}}}{0.03}\right)}^{-1}.$$(9)

Fig. 1 Width of the outflow region as a function of the planetary mass (Eq. (12)).

2.5.2 Definition of the width of the outflow region

The dust velocity is significantly perturbed at the edges of the outflow region. Following Kuwahara & Kurokawa (2024), we refer to the outflow region as the region where the x-component of the gas velocity is dominantly perturbed. The dimensionless x-coordinate of the edge of the outflow region is given by (Kuwahara & Kurokawa 2024) $$w_{\text{out}}^{\pm }=\pm \min \left(\frac{2}{3}\left(1\mp {ℳ}_{\text{hw}}\right),w_{\text{HS}}+\frac{2}{3}{ℳ}_{\text{hw}}\right),$$(10)

where $$w_{\text{HS}}=\frac{1.05\sqrt{m}+3.4m^{7/3}}{1+2m^2}$$(11)

is the half-width of the horseshoe region (Jiménez & Masset 2017). From Eq. (10), the width of the outflow region is given by $${\mathcal{W}}_{\text{out\hspace{0.17em}}}\equiv w_{\text{out\hspace{0.17em}}}^{+}-w_{\text{out\hspace{0.17em}}}^{-}=\min \left(2w_{\text{HS}}+\frac{4}{3}{ℳ}_{\text{hw}},\frac{4}{3}\right).$$(12)

We plotted Eq. (12) in Fig. 1. The width of the outflow region increases with the planetary mass when m ≲ 0.3 and converges at m ≳ 0.3.

2.5.3 Definition of the width of the dust ring and gap

Following Paper I, hereafter we refer to the regions where dust is depleted and accumulated as “dust gap” and “dust ring”, respectively. We numerically calculated the dust gap width by $${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)=x_{\text{gap,out\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)-x_{\text{gap,in\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t),$$(13)

where the superscript “num” represents the value obtained from numerical simulations. In the above equation, $x_{\text{gap,in\hspace{0.17em}}}^{\text{num}}(t)$ and $x_{\text{gap,out}}^{\text{num}}(t)$ are the edges of the dust gap where Σ_d(x, t) reaches the following value: (Dong & Fung 2017)² $${\text{Σ}}_{\text{d}}\left(x_{\text{gap\hspace{0.17em}},\text{\hspace{0.17em}in\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)\right)={\text{Σ}}_{\text{d}}\left(x_{\text{gap\hspace{0.17em}},\text{\hspace{0.17em}out\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)\right)=\sqrt{{\text{Σ}}_{\text{d},0}\times \overline{{\text{Σ}}_{\text{d},\min }}(t)}.$$(14)

We defined $\overline{{\text{Σ}}_{\text{d},\min }}(t)$ as the average of the minimum dust surface density inside and outside the planetary orbit: $$\overline{{\text{Σ}}_{\text{d},\min }}(t)\equiv \frac{{{\text{Σ}}_{\text{d},\min }(t)|}_{x<0}+{{\text{Σ}}_{\text{d},\min }(t)|}_{x\ge 0}}{2}.$$(15)

We numerically calculated the dust ring width as $${\mathcal{W}}_{\text{ring\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)=x_{\text{ring,out\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)-x_{\text{ring,in\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t).$$(16)

In Eq. (16), $x_{\text{ring,\hspace{0.17em}in\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)\text{\hspace{0.17em}and\hspace{0.17em}}{x}_{\text{ring,out\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)$ are the edges of the dust ring where Σ_d(x, t) reaches the geometric mean of the maximum and initial values: $${\text{Σ}}_{\text{d}}\left(x_{\text{ring,in\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)\right)={\text{Σ}}_{\text{d}}\left(x_{\text{ring,out\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)\right)=\sqrt{{\text{Σ}}_{\text{d},0}\times {{\text{Σ}}_{\text{d},\max }|}_{x\ge 0}(t)}.$$(17)

For the definition of the dust ring, we only focus on the dust accumulation outside the planetary orbit. The definitions of the widths of the dust ring and gap were plotted in Fig. 2.

Fig. 2 Definition of the widths of the dust ring and gap. The red and yellow shaded regions show the numerically calculated widths of the dust ring and gap. The circles mark Σd,max and Σd,min.

2.5.4 Definition of the depth of the dust gap

We defined the dust gap depth as the contrast between the minimum and maximum dust surface densities (Fig. 2; Huang et al. 2018; Zhang et al. 2018). We numerically calculated the dust gap depth as $${\delta }_{\text{gap\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)\equiv \frac{{\text{Σ}}_{\text{d},\min }(t)}{{\text{Σ}}_{\text{d},\max }(t)}.$$(18)

It should be noted that Eq. (18) may represent the amplitude of the dust ring if turbulent diffusion smooths the dust gap and then only the dust ring forms. Although it can differ from the intuitive definition of the dust gap depth, we consistently use Eq. (18) for measuring the dust gap depth.

Fig. 3 Perturbation of the planet-induced gas flow on the radial velocity of dust and the dust surface density. We set m = 0.1, ℳhw = 0.03, St = 10−3, and αdiff = 10−4. Top: gas flow structure at the meridian plane. The color contour represents the gas velocity in the x-direction averaged in the y-direction within the calculation domain of hydrodynamical simulation, 〈vx,g〉y. The vertical dotted lines represent the x-coordinates of the edges of the outflow region, ${\omega }_{\text{out\hspace{0.17em}}}^{\pm }$. Middle: perturbed radial drift velocity of dust. The horizontal dashed line represents vdrift. Bottom: time evolution of the dust surface density. The gray dashed line corresponds to the steady-state dust surface density. The circle and triangle symbols denote the location of the numerically calculated edges of the dust gap and ring.

3 Numerical results: time evolution of Σ_d(t)

In this section we present the numerical results of this study, comparing them with the steady-state solutions Paper (I). We describe the physical processes of the time evolution of the dust surface density. Based on the results presented in this section, we introduce the semi-analytic models of the widths of the dust ring and gap and the depth of the dust gap in Sect. 4.4.

We describe the behavior of the time evolution of Σ_d(t) influenced by the planet-induced gas flow. Given the wide parameter spaces in our study, we first show the numerical results for a fiducial parameter set with α_diff = 10⁻⁴, m = 0.1, ℳ_hw = 0.03, and St = 10⁻³ (Sect. 3.1). We then show the dependence of Σ_d(t) on the turbulent parameter (Sect. 3.2), the planetary mass (Sect. 3.3), the Mach number of the headwind (Sect. 3.4), and the Stokes number (Sect. 3.5), respectively.

3.1 Fiducial case

The outflow of the gas at the midplane induced by low-mass planets perturbs the radial drift velocity of dust, causing the dust rings and gaps (Fig. 3). Figure 3 shows the perturbations of low- mass planets on the gas and dust, where we set α_diff = 10⁻⁴, m = 0.1, ℳ_hw = 0.03, and St = 10⁻³, as the fiducial parameter set. In Fig. 3, the top, the middle, and the bottom panels show the velocity field of the gas at the meridian plane, the radial drift velocity of dust influenced by the gas flow, and the dust surface density, respectively.

The dust surface density, which has initially a flat profile, changes with time, and then reaches a steady state within t ≲ 10⁶ (Fig. 3). The outflow of the gas at the midplane perturbs the radial drift velocity of dust, causing positive and negative peaks in the profile of 〈v_d〉 (the middle panel of Fig. 3). The radially- outward (inward) outflow of the gas inhibits (enhances) the radial drift of dust. The dust surface density decreases with time around the planetary orbit creating a dust gap, while it increases with time outside the planetary orbit creating a dust ring. The dust surface density only changes by a factor of ∼2 in Fig. 3 due to efficient dust diffusion. Given a characteristic length of a perturbation is set by ℒ = 𝒲_out, the drift and diffusion timescales are given by t_drift ~ 1.2 × 10⁴ and t_drift ~ 5.2 × 10³, resulting in a diffusion-dominated regime t_diff < t_drift.

The locations of the edges of the dust gap are determined by those of the edges of the outflow region, $x={\omega }_{\text{out\hspace{0.17em}}}^{\pm }$ (Eq. (10); the vertical dotted lines in Fig. 3). Thus, the dust gap widths can be estimated by 𝒲_out (Eq. (12)). The location of the inner edge of the dust ring is set by $x={\omega }_{\text{out\hspace{0.17em}}}^{\pm }$ and hardly changes with time. The outer edge of the dust ring moves with time due to diffusion (Fig. 3).

Fig. 4 Time evolution of the dust surface density in low-turbulence disks. We set m = 0.1, ℳhw = 0.03, St = 10−3, and αdiff = 10−5. The horizontal axis is on a log scale, whose range is extended to x ∈ [−100, 5]. The vertical dashed lines show the analytic model of the location of the inner edge of the dust gap, which moves with vdrift (Sect. 4.4).

3.2 Dependence of Σ_d (t) on the turbulent parameter

A perturbation to Σ_d (t) due to the planet-induced gas flow strongly depends on the turbulent parameter. Figure 4 shows the time-dependent Σ_d(t) in the low-turbulence disk (α_diff = 10⁻⁵). The planetary mass, the Mach number of the headwind, and the Stokes number are the same as in the fiducial case. Compared to the fiducial case, Σ_d (t) changes more significantly, because a steep gradient of Σ_d needs to achieve a steady state due to inefficient dust diffusion (t_diff > t_drift in Fig. 4). In Fig. 4, Σ_d(t) increases (decreases) by ∼3–4 orders of magnitude. The dust surface density reaches the steady state after t > 10⁸.

During the time evolution Σ_d(t) changes in a complex manner. The profile of Σ_d(t) deviates significantly from the steadystate solution. The dust accumulates over time at $x\gtrsim {\omega }_{\text{out\hspace{0.17em}}}^{+},$, which is similar to the fiducial case. At $x\lesssim {\omega }_{\text{out\hspace{0.17em}}}^{+}$, the dust surface density drops significantly in the early stage of the time evolution, t ≲ 10⁴–10⁵. Due to inefficient dust diffusion, the dust hardly leak out to the inside of the planetary orbit. As a result, at the early stage of the time evolution, the minimum value of Σ_d(t) can be orders of magnitude smaller than that of the steady-state solution.

Moreover, we found that the dust gap expands with time, and then the dust is depleted in a wide range of $x\lesssim {\omega }_{\text{out\hspace{0.17em}}}^{+}$ when m ≳ 0.1 (Fig. 4; see also Fig. C.1 for different planetary masses; available on Zenodo³). We found that the inner edge of the dust gap moves with |v_drift|. In this case, the dust gap width can no longer estimated by 𝒲_out (Eq. (12)). Because we assumed a constant supply of dust from the outer region of the disk, the dust slowly leaks to the inside of the planetary orbit, and then Σ_d (t) at $x\lesssim {\omega }_{\text{out\hspace{0.17em}}}^{+}$ increased in the late stage of the time evolution, t ≳ 10⁸.

For illustrative purposes, we also display these time sequences in a two-dimensional (2D) plane in Fig. 5, which were generated from the results of 1D calculations shown in Fig. 4. The dust cavity, the gap as wide as the planet’s orbital distance, forms after t ≥ 10⁶ . The cavity-opening timescale can be estimated by Eq. (7) with ℒ = r_p. We note that our model for the dust surface density evolution ignores the effect of curvature (Eq. (5)).

A disk with a dust cavity formed by the gas-flow mechanism could be observed as a transition disk (Francis & van der Marel 2020), although the dependence of the cavity evolution time on m, ℳ_hw, St, and α_diff remains unclear. We speculate that the dust cavity is filled on long timescale (≳10⁶–10⁷ , corresponding to > 5–50 Myr at 10 au; Fig. 6). Since the cavity is filled by diffusion of dust at the ring, we considered that the cavity-filling timescale would be proportional to the formation timescale of the dust ring (Eq. (35) of Paper I): τ_cav,fill ~ M_ring/Ṁ_dust, where M_ring is the mass of the steady-state dust ring and Ṁ_dust is the radial inward mass flux of dust. We discuss the implications for transition disks in Sect. 5.3.1.

Fig. 5 Time evolution of Σd(t) in the two-dimensional plane. We set m = 0.1, ℳhw = 0.03, St = 10−3, and αdiff = 10−5. These images were generated based on the results of 1D calculations assuming an axisymmetric dust distribution, neglecting disk curvature. The axes are normalized by the planet location, rp, calculated by X = Y = (r − rp)/hp, where r is the radial coordinate centered at the star and hp is the disk aspect ratio at r = rp. We set hp = 0.05.

3.3 Dependence of Σ_d (t) on the planetary mass

When higher-mass planets are assumed, we find deeper dust gaps and higher concentrations of dust in a ring. In Fig. 7, the Mach number of the headwind, the Stokes number, and the turbulent parameter are the same as in the fiducial case: ℳ_hw = 0.03, St = 10⁻³, and α_diff = 10⁻⁴. The location of the outer edge of the dust gap, $x_{\text{gap,out\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}~{\omega }_{\text{out\hspace{0.17em}}}^{+}$, hardly changes when m ≳ 0.3, at which point ${\omega }_{\text{out\hspace{0.17em}}}^{+}$ converges (Eq. (10)). The dust gap depths converge at m ≳ 0.3. This is because the outflow speed in the x-direction has a peak at m ~ 0.3 and then the influence of the gas flow on the dust motion saturates (Kuwahara & Kurokawa 2024). In Paper I, where we only considered m = 0.03, 0.1, and 0.3, we set the condition of the dust ring and gap formation as m ≳ 0.1. However, we found that even low-mass planets with m = 0.05 can generate dust rings and (or) gaps (zoom-in panels of Fig. 7, see also Fig. C.1 for α_diff = 10⁻⁵).

Fig. 6 Cavity-filling timescale. We only focus on m ≥ 0.1 in which the dust cavities form.

3.4 Dependence of Σ_d (t) on the Mach number of the headwind

The amplitude of the dust ring becomes higher when the larger ℳ_hw is assumed (Fig. 8). The flow speed of the radially-outward outflow increases with ℳ_hw (Kuwahara & Kurokawa 2024), which leads to higher concentrations of dust into a ring outside the planetary orbit.

3.5 Dependence of Σ_d (t) on the Stokes number

In this section we focus on the dependence of Σ_d(t) on the Stokes number. In Fig. 9, the planetary mass, the Mach number of the headwind, and the turbulent parameter are the same as in the fiducial case: m = 0.1, ℳ_hw = 0.03, and α_diff = 10⁻⁴.

Paper I found that the deeper dust gaps and higher concentrations of dust into a ring can be seen for the smaller Stokes numbers in a steady state, because smaller dust particles are more sensitive to the gas flow. This is successfully reproduced in Fig. 9d. The time required to reach the steady state is shorter when the larger Stokes number is assumed (Figs. 9a–c, see also Fig. C.2 for α_diff = 10⁻⁵), because the drift timescale of dust is shorter for the larger Stokes numbers (t_drift ∝ St⁻¹ ; Eq. (7)).

Fig. 7 Dependence of Σd(t) on the planetary mass. We set ℳhw = 0.03, St = 10−3, and αdiff = 10−5. The vertical dotted lines correspond to |x| = 4/3 (the x-coordinate of the edge of the outflow region for m ≳ 0.3; Eq. (12)). The small panels above the upper left corners of panels a-c are the zoomed-in views for m = 0.03, 0.05, and 0.07.

Fig. 8 Dependence of Σd(t) on the Mach number of the headwind. We set m = 0.2, St = 10−3, and αdiff = 10−4. We varied the Mach number of the headwind in each panel, ℳhw = 0.01 (top) and ℳhw = 0.1 (bottom).

3.6 Summary of the parameter dependence

We summarize the dependence of Σ_d(t) on α_diff, m and St for a fixed ℳ_hw and the time in Fig. 10. We set ℳ_hw = 0.03 and t = 10⁵ in Fig. 10. A perturbation to Σ_d(t) is stronger when the smaller α_diff, the higher-mass planets, or the smaller St are assumed.

Figure 11 is a contour plot of the dust gap depth as a function of the planetary mass and the Stokes number for a fixed ℳ_hw and time (ℳ_hw = 0.03 and t = 10⁵), showing that the dust gap forms when m ≳ 0.05 and the dust gap deepens as the planetary mass increases. At t = 10⁵, the dust gap depth has a peak at St = 10⁻³.

4 Properties of dust rings and gaps

This section shows the widths of the ring and gap and depth of the gap. We first show the numerical results in Sects. 4.1–4.3. We then introduce semi-analytic models of dust rings and gaps in Sect. 4.4 based on the obtained numerical results.

4.1 Numerically calculated dust gap width

We found that the dust gap width either stays constant or expands with time (Fig. 12). In general, once the dust gaps form, their widths do not change significantly with time in the moderateturbulence disks (α_diff = 10⁻⁴; Fig. 12a), because Σ_d(t) decreases only within the outflow region (Fig. 3). The dust gap widths increase with the planetary mass when m ≲ 0.3, and converge at m ≳ 0.3. These numerical results can be explained by the dependence of 𝒲_out on the planetary mass, which is independent of time. It increases with the planetary mass when m ≲ 0.3, and has a constant value at m ≳ 0.3 (Fig. 1).

The dust gap keeps expanding with time when m ≳ 0.1 after t ≳ 10⁴ in the low-turbulence disk (α_diff = 10⁻⁵; Fig. 12b). Since the inner edge of the dust gap moves with |v_drift| (Fig. 4), the width of the expanding dust gap is independent of the planetary mass. We note that the semi-analytic model of the dust gap width (the dashed lines in Fig. 12) will be introduced in Sect. 4.4.

4.2 Numerically calculated dust gap depth

The dust gaps deepen initially with time, and then their depths ${\delta }_{\text{gap\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)$ converge after t ≳ 10³–10⁴ (Fig. 13). Initially, ${\delta }_{\text{gap\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}(t)$ decreases because the dust surface density decreases at $x\lesssim {\omega }_{\text{out\hspace{0.17em}}}^{+}.$. As Σ_d(t) stops decreasing at $x\lesssim {\omega }_{\text{out\hspace{0.17em}}}^{+}$ or a decrease in Σ_d(t) at $x\lesssim {\omega }_{\text{out\hspace{0.17em}}}^{+}$ balances an increase in Σ_d(t) at $x\gtrsim {\omega }_{\text{out\hspace{0.17em}}}^{+}$ after t ≳ 10³–10⁴, the dust gap depth eventually becomes constant (Fig. 4). The dust gaps deepen with the planetary mass when m ≲ 0.3 and their depths converge at m ≳ 0.3 (Fig. 14), because the outflow speed has a peak at m ∼ 0.3 and, consequently, the perturbation of the gas flow on the dust motion saturates (Kuwahara & Kurokawa 2024). We note that the semi-analytic model of the dust gap depth (the dashed lines in Fig. 13) will be introduced in Sect. 4.4.

Fig. 9 Dependence of Σd(t) on the Stokes number. We set m = 0.1, ℳhw = 0.03, and αdiff = 10–5. The vertical dotted lines correspond to $x={\omega }_{\text{out\hspace{0.17em}}}^{\pm }$.

Fig. 10 Summary of the parameter dependence of Σd(t). We set ℳhw = 0.03 and t = 105. We varied the planetary mass, the Stokes number, and the turbulent parameter.

4.3 Numerically calculated dust ring width

The dust ring widths increase with time due to diffusion and then reach a steady state (Fig. 15). The dust ring widths have the radial extent of ≲1–10 times the gas scale height, which is weakly dependent on the planetary mass. Figure 15 summarizes the parameter dependence of the dust ring width at a certain time, showing the following trends. (1) The dust ring width decreases as St increases (Figs. 15a and b). (2) The dust ring width increases as α_diff increases (Figs. 15b and c). (3) The dust ring width decreases as ℳ_hw increases (Figs. 15c and d). These trends suggest that the dust ring width is proportional to the length where the drift timescale coincides with the diffusion timescale: ${ℒ}_{\text{eq}}\propto {\alpha }_{\text{diff}}{\text{St}}^{-1}{ℳ}_{\text{hw}}^{-1}$. We note that the semi-analytic model of the dust ring width (the dashed lines in Fig. 15) will be introduced in Sect. 4.4.

Fig. 11 Contour plot of the dust gap depth as a function of the planetary mass and the Stokes number. We set ℳhw = 0.03 and t = 105. We varied the turbulent parameter in each panel, αdiff = 10–4 (panel a) and αdiff = 10–5 (panel b).

4.4 Semi-analytic models of dust rings and gaps

Based on the obtained numerical results in Sect. 3, we introduce semi-analytic models of the width of the dust ring and gap and the depth of the dust gap as functions of m, ℳ_hw, St, α_diff, and t. Since a significant perturbation to Σ_d(t) due to the planet- induced gas flow appears only when the smaller Stokes numbers were assumed, we restrict our attention to the limited range of the Stokes number, St ≲ 10^–3.

Section 4.4.1 introduces the semi-analytic model of the dust gap width. By fitting of the numerical results, we derived a criterion for the dust gap width which distinguish between the temporally constant and expanding gaps, and then described the dust gap widths in each case. In Sect. 4.4.2, we considered that the time evolution of the dust gap depth is described by a logistic differential equation. Using an analytical solution to the logistic equation combined with the fitting of numerical results, we obtained the semi-analytic model of the dust gap depth. In Sect. 4.4.3, we considered the time evolution of the dust ring width is described by a sigmoid curve with a steady-state dust ring width as an asymptote. By fitting the sigmoid curve to the numerical results, we obtained the semi-analytic model of the dust ring width.

Fig. 12 Time evolution of the dust gap width for different planetary masses. We fixed the Stokes number St = 10–3 and the Mach number ℳhw = 0.03, and set αdiff = 10–4 in panel a and αdiff = 10–5 in panel b. The solid lines with the circle symbols and the dashed lines are the numerically calculated and the semi-analytic dust gap widths, respectively (Eq. (26); Sect. 4.4). We note that in panels a and b, the numerically calculated dust gap width for m = 0.03 is not shown because we obtained ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{num\hspace{0.17em}}}=0$.

4.4.1 Dust gap width

As mentioned in Sect. 4.1, the dust gap is either constant or expanding with time. The temporally constant dust gap width can be modeled by the width of the outflow region, W_out (Eq. (12)). When the dust gap expands with time, the inner edge of the dust gap is set by $x_{\text{gap,\hspace{0.17em}in\hspace{0.17em}}}(t)=-\left|{\upsilon }_{\text{drift\hspace{0.17em}}}\right|t\text{\hspace{0.17em}or\hspace{0.17em}}{\omega }_{\text{out\hspace{0.17em}}}^{-}$, whichever smaller (Fig. 4).

Considering the dust surface density within the outflow region, we construct a semi-analytic model for the dust gap width, ${\mathcal{W}}_{\text{gap}}^{\text{SA}}(t)$. We expect that the dust gap width keeps constant, ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{SA}}(t)={\mathcal{W}}_{\text{out}}$, when the diffusive flux of dust dominates the time evolution of the dust surface density, while we assume that the dust gap expands with time when the advective flux dominates. We determine the diffusion-dominated or advection- dominated regime by comparing Σ_d(x, t) at the gap location with a certain critical value, Σ_crit. Given the balance between the advective and diffusive flux of dust, we derive the critical dust surface density: $$\langle v_{\text{d}}\rangle {\text{Σ}}_{\text{d}}=\mathcal{D}\frac{\partial {\text{Σ}}_{\text{d}}}{\partial x}.$$(19)

Here we focus on the limited region, ${\omega }_{\text{out\hspace{0.17em}}}^{-}<x<{\omega }_{\text{out\hspace{0.17em}}}^{+}$, in which the radial drift velocity of dust is approximately given by ⟨v_d⟩ ∼ v_dfrit (the middle panel of Fig. 3). We set ⟨v_d⟩ = v_drift for simplicity in Eq. (19). We then integrate Eq. (19) over a range of ${\mathcal{W}}_{\text{out\hspace{0.17em}}}={\omega }_{\text{out\hspace{0.17em}}}^{+}-{\omega }_{\text{out\hspace{0.17em}}}^{-}$ (Eq. (12)) and obtain $$\ln \left(\frac{{\text{Σ}}_{\text{d}}\left(w_{\text{out\hspace{0.17em}}}^{+}\right)}{{\text{Σ}}_{\text{d}}\left(w_{\text{out\hspace{0.17em}}}^{-}\right)}\right)=\frac{v_{\text{drift\hspace{0.17em}}}}{\mathcal{D}}{\displaystyle {\int }_{w_{\text{out\hspace{0.17em}}}^{\text{-}}}^{w_{\text{out\hspace{0.17em}}}^{+}}\text{d}}x.$$(20)

Equation (20) gives $${\text{Σ}}_{\text{d}}\left(w_{\text{out\hspace{0.17em}}}^{+}\right)={\text{Σ}}_{\text{d},0}\exp \left[-\frac{2\text{St}{ℳ}_{{}_{\text{hw}}}{\mathcal{W}}_{{}_{\text{out}}}}{{\alpha }_{\text{diff\hspace{0.17em}}}}\right]\equiv {\text{Σ}}_{\text{crit\hspace{0.17em}}},$$(21)

where we set ${\text{Σ}}_{\text{d}}\left({\omega }_{\text{out\hspace{0.17em}}}^{-}\right)={\text{Σ}}_{\text{d},0}$ and assume 1 + St² ≃ 1. The diffusive (advective) flux of dust dominates the time evolution of the dust gap when Σ_d(x, t) > ∑_crit (Σ_d(x, t) < ∑_crit) within the limited region, ${\omega }_{\text{out\hspace{0.17em}}}^{-}<x<{\omega }_{\text{out\hspace{0.17em}}}^{+}$. We compared the time evolution of the dust surface density with Σ_crit in Fig. A.2.

Practically, we consider that the dust gap stays constant when ${\text{Σ}}_{\min }^{\text{global\hspace{0.17em}}}\ge {\text{Σ}}_{\text{crit\hspace{0.17em}}}$ and expands with time when ${\text{Σ}}_{\min }^{\text{global\hspace{0.17em}}}\ge {\text{Σ}}_{\text{crit\hspace{0.17em}}}$, where ${\text{Σ}}_{\min }^{\text{global\hspace{0.17em}}}$ is the global minimum of the time-dependent dust surface density at the gap location: $${\text{Σ}}_{\min }^{\text{global\hspace{0.17em}}}\equiv \underset{t>0}{\min }{\text{Σ}}_{\text{d},\min }(t).$$(22)

We find that ${\text{Σ}}_{\min }^{\text{global\hspace{0.17em}}}$ can be fitted by (Appendix B) $${\text{Σ}}_{\min }^{\text{global\hspace{0.17em}}}=\min \left(1,{\text{Σ}}_{\min }^{\text{fit}}\right),$$(23)

where $${\text{Σ}}_{\min }^{\text{fit}}={10}^{{\mathcal{S}}_{\min }^{\text{fit}}\left({\alpha }_{\text{diff\hspace{0.17em}}},m\right)},$$(24) $${\mathcal{S}}_{\min }^{\text{fit}}\left({\alpha }_{\text{diff\hspace{0.17em}}},m\right)=-0.37{\left(\frac{{\alpha }_{\text{diff\hspace{0.17em}}}}{{10}^{-4}}\right)}^{-1.1}\times \text{erf}\left(3.2\times {10}^2{\left(\frac{{\alpha }_{\text{diff\hspace{0.17em}}}}{{10}^{-4}}\right)}^{-0.17}{m}^{2.8}\right),$$(25)

with erf being the error function (erf(m) → 1 when m → ∞).

In summary, the semi-analytic formula of the dust gap width is given by $${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{SA\hspace{0.17em}}}(t)=\{\begin{array}{l}{\mathcal{W}}_{\text{out\hspace{0.17em}}}\left({\text{Σ}}_{\min }^{\text{global\hspace{0.17em}}}\ge {\text{Σ}}_{\text{crit\hspace{0.17em}}}\right),\hfill \\ \max \left({\mathcal{W}}_{\text{out\hspace{0.17em}}},w_{\text{out\hspace{0.17em}}}^{+}-\left|v_{\text{drift\hspace{0.17em}}}\right|t\right)\left({\text{Σ}}_{\min }^{\text{global\hspace{0.17em}}}<{\text{Σ}}_{\text{crit\hspace{0.17em}}}\right).\hfill \end{array}$$(26)

We plotted Eq. (26) in Fig. 12 with the dashed line (see also Fig. C.3). In Fig. 12, Eq. (26) predicts that the dust gaps keep expanding with time when m ≳ 0.1–0.2 at t ≳ 10⁴, which is consistent with the numerical results in the low-turbulence disk (α_diff = 10^–5; Fig. 12b). When α_diff = 10^–4, Eq. (26) fails to reproduce the numerical results for m ≥ 0.2 which are constant with time (Fig. 12a). We speculate that this deviation is caused by the assumption of ⟨v_d⟩ ∼ v_dfrit in Eq. (19). Nevertheless, we use Eq. (19) as it reproduces an overall trend in the planetary-mass dependence.

Fig. 13 Time evolution of the dust gap depth for different planetary masses. We fixed the Stokes number St = 10–3 and the Mach number ℳhw = 0.03, and set αdiff = 10–4 in panel a and αdiff = 10–5 in panel b. The solid lines with the circle symbols and the dashed lines are the numerically calculated and the semi-analytic dust gap depths, respectively (Eq. (28); Sect. 4.4).

Fig. 14 Dust gap depth as a function of planetary mass at t = 106. We fixed the Stokes number St = 10–3 and set αdiff = 10–4 in panel a and αdiff = 10–5 in panel b. The solid lines with the circle symbols and the dashed lines are the numerically calculated and the semi-analytic dust gap depths, respectively (Eq. (28); Sect. 4.4).

Fig. 15 Time evolution of the dust ring width for different planetary masses. The assumed parameters (ℳhw, St, and αdiff) are shown at the top of each panel. The solid lines with the circle symbols and the dashed lines are the numerically calculated and the semi-analytic dust ring widths, respectively (Eq. (33); Sect. 4.4).

4.4.2 Dust gap depth

As described in Sect. 4.2, the dust gap depth deepens with time and has a lower limit. We assume that the time evolution of the dust gap depth obeys the following equation⁴: $${\delta }_{\text{gap\hspace{0.17em}}}^{SA}(t)={\delta }_{\infty }{\left[1-\left(1-\frac{{\delta }_{\infty }}{{\delta }_0}\right)e^{-t/{\tau }_{\text{gap\hspace{0.17em}}}}\right]}^{-1},$$(28)

where ${\delta }_{\text{gap\hspace{0.17em}}}^{\text{SA}}(t)$ is the semi-analytic model of the dust gap depth, δ_∞ is the steady-state dust gap depth, ${\delta }_0={\delta }_{\text{gap\hspace{0.17em}}}^{\text{SA}}(0)\equiv 1$, and τ_gap is the characteristic timescale. Equation (28) shows that ${\delta }_{\text{gap\hspace{0.17em}}}^{\text{SA}}(t)$ decreases with time and then approaches the steady-state value, ${\delta }_{\text{gap\hspace{0.17em}}}^{\text{SA}}(t)\to {\delta }_{\infty }$, at t >> τ_gap. We define $${\tau }_{\text{gap}}\equiv \min \left(t_{\text{drift}},t_{\text{diff}}\right),$$(29)

where we set ℒ = 0.41 × W_out in both t_drift and t_diff (Eqs. (7) and (8)). The coefficient of 0.41 is determined by the least-squares fitting of numerical results. The characteristic timescale τ_gap is a function of m, ℳ_hw, St, and α_diff, having on the order of ~ 10³ –10⁴ in our parameter sets.

We find that the steady-state dust gap depth can be fitted by (Appendix B) $${\delta }_{\infty }=\min \left(1,{\delta }_{\infty }^{\text{fit}}\right),$$(30)

where $${\delta }_{\infty }^{\text{fit}}={10}^{{\mathcal{S}}_{\infty }^{\text{fit}}\left({\alpha }_{\text{dif}},m\right)},$$(31) $${\mathcal{S}}_{\infty }^{\text{fit}}\left({\alpha }_{\text{diff}},m\right)=-0.63{\left(\frac{{\alpha }_{\text{diff}}}{{10}^{-4}}\right)}^{-1.1}\times \mathrm{erf}\left(4.2\times {10}^2{\left(\frac{{\alpha }_{\text{diff}}}{{10}^{-4}}\right)}^{0.022}{m}^{2.8}\right).$$(32)

We plotted Eq. (28) in Figs. 13 and 14 (see also Fig. C.4) with the dashed line. Although Eq. (28) does not completely reproduce the numerical results, it shows good agreement with the numerical result, in particular when t ≳ τ_gap ∼ 10³–10⁴. When t ≲ τ_gap , Eq. (28) only agrees with the numerical results of m < 0.1. We speculate that the deviation is caused by the assumption of τ_gap , which is set by the drift or the diffusion timescale with ℒ ∝ W_out, whichever smaller (Eq. (29)). However, the radial drift speed of dust deviates from the unperturbed value within the outflow region, which changes t_drift.

4.4.3 Dust ring width

As described in Sect. 4.3, the dust ring width increases with time and then reaches a steady state. We assume that the time evolution of the dust ring width is described by the following sigmoid curve: $${\mathcal{W}}_{\text{ring\hspace{0.17em}}}^{\text{SA}}(t)={\mathcal{W}}_{\text{ring,\hspace{0.17em}},\infty }^{\text{fit}}\left(1-\frac{1}{1+{\left(t/{\tau }_{\text{ring\hspace{0.17em}}}\right)}^q}\right).$$(33)

Here ${\mathcal{W}}_{\text{ring\hspace{0.17em}}}^{\text{SA}}(t)$ is the semi-analytic model of the dust ring width, ${\mathcal{W}}_{\text{ring,}\infty }^{\text{fit\hspace{0.17em}}}$ is the fitting formula for the steady-state dust ring width, τ_ring is the characteristic timescale, and q = 0.42 (Appendix B). Numerical results showed that ${\mathcal{W}}_{\text{ring,\hspace{0.17em}},\infty }^{\text{fit}}$ would be proportional to L_eq (Eq. (9)), and, consequently, α_diff (Sect. 4.3), but we find that the dependence is weaker than ${\mathcal{W}}_{\text{ring,\hspace{0.17em}}\infty }^{\text{fit\hspace{0.17em}}}\propto {\alpha }_{\text{diff\hspace{0.17em}}}$. We find that ${\mathcal{W}}_{\text{ring,\hspace{0.17em}}\infty }^{\text{fit\hspace{0.17em}}}$ can be fitted by (Appendix B): $${\mathcal{W}}_{\text{ring,\hspace{0.17em}}\infty }^{\text{fit}}=0.63{\left(\frac{{\alpha }_{\text{diff}}}{{10}^{-4}}\right)}^{-0.65}\times {ℒ}_{\text{eq}}.$$(34)

The dust rings expand due to dust diffusion. Thus, we define $${\tau }_{\text{ring\hspace{0.17em}}}\equiv \frac{{\left({\mathcal{W}}_{\text{ring,\hspace{0.17em}}\infty }^{\text{fit}}\right)}^2}{{\alpha }_{\text{diff\hspace{0.17em}}}}.$$(35)

We plotted Eq. (33) in Fig. 15 (see also Fig. C.5) with the dotted lines, which show good agreement with the numerical results.

Fig. 16 Same as Fig. 14, but for St = 10–2.

4.4.4 Caveat

So far we have considered the regime in which the dust is tightly coupled with the gas, St ≲ 10^–3. Since we developed our semi-analytic models by fitting the numerical results with St = 10^–3 (Appendix B), our semi-analytic models would be invalid when St ≳ 10^–2. Figure 16 compares our semi-analytic model for the dust gap depth with the numerical result when St = 10^–2, showing a significant deviation appears in particular when α_diff = 10^–5.

Fig. 17 Time-dependent dust surface density with particle sizes of s = 3.7 mm perturbed by an Earth-like planet (Mp = 0.66M⊕) at 1 au. Numerical simulations were conducted in dimensionless units. We set m = 0.1, ℳhw = 0.03, St = 10–3, and varied the turbulent parameter in each panel, αdiff = 10–5 (panel a) and αdiff = 10–4 (panel b). For the assumed parameter set, the pebble isolation masses are given by Miso = 2.8 M⊗ (panel a) and 3 M⊗ (panel b; Eq. (A.14)).

5 Discussion

5.1 Time evolution of Σ_d (t) with dimensional units

To facilitate the interpretation of our results, we rewrite our numerical results with dimensional units. Assuming that a typical steady accretion disk model (Ida et al. 2016), we convert the dimensionless quantities into the dimensional ones. Appendix A describes the method for the conversion. We set the orbital radius of the planet as r_p = 1 au or 50 au, at which the Mach number of the headwind has the value of ℳ_hw = 0.03 or 0.1 (Eq. (A.10)).

Here we show the time-dependent dust surface density Σ_d(t) with sizes of s ≈ 4 mm-sized particles (St = 10^–3) perturbed by gas flow induced by an Earth-like planet at 1 au (Sect. 5.1.1). We also show solids with sizes of s ≈ 0.2 mm (St = 3 × 10^–3) perturbed by a Neptune-like planet at 50 au (Sect. 5.1.2). The dust size was chosen to be consistent with nonsticky slicate grains inside the H₂O snow (≲ a few au) and with nonsticky icy CO₂- covered grains outside the CO₂ snow line located approximately outside 10 au (Musiolik et al. 2016a,b). This limits particle sizes to ~2 mm at ~1 au and ~0.4 mm at 50 au (Okuzumi & Tazaki 2019).

5.1.1 Earth-like planet at 1 au

Figure 17 shows the evolution of the solid surface density Σ_d(t) of ~4 mm sized particles around an Earth-like planet with a mass of ~0.7 M_⊗ embedded at 1 au in a disk with low midplane turbulence (α_diff = 10^–5; Fig. 17a) and in a disk with moderater midplane turbulence (α_diff = 10^–4; Fig. 17b). In the low-turbulence disk, the dust is depleted by ~2 orders of magnitude within 1 Myr at <1 au (Fig. 17b). A significant amount of dust concentrates into a very narrow ring whose width is less than 0.1 au at the planet location. In the narrow ring, Σ_d (t) increases by ~10² times the initial value. In contrast, in the moderate-turbulence disk (α_diff = 10^–4; Fig. 17a), the Earthmass planets do not create a significant dust depletion nor concentration for the assumed parameters. A shallow dust gap appears within 0.1 Myr at <1 au, but it is smoothed within 1 Myr. Only a narrow dust ring whose width is ~ 0.1 au remains outside the planet location at 1 Myr.

The assumed planetary mass (M_p ≃ 0.7 M_⊗) falls below the pebble isolation mass (M_iso ≈ 3 M_⊗; Lambrechts et al. 2014; Bitsch et al. 2018) at which a growing planet opens a shallow gas gap and then the pebble flux is highly reduced inside the planetary orbit. Even planets with masses below the pebble isolation mass can affect significantly Σ_d(t). When the planetary mass exceeds m ≳ 0.1 (M_p ≳ 0.7 M_⊗ at 1 au), the subsequent growth of the planet would be suppressed and planets remain small in the terrestrial planet forming-region.

Fig. 18 Time-dependent dust surface density with particle sizes of s = 0.23 mm for a Neptune-like planet (Mp = 13 M⊗) at 50 au. The numerical simulations were conducted in dimensionless units. We set m = 0.1, ℳhw = 0.1, St = 3 × 10–3, and varied the turbulent parameter in each panel, αdiff = 10–5 (panel a) and αdiff = 10–4 (panel b). For the assumed parameter set, the pebble isolation masses are given by Miso = 56 M⊗ (panel a) and 60 Mæ (panel b; Eq. (A.14)).

5.1.2 Neptune-like planet at 50 au

Figure 18 shows a situation where a Neptune-like planet, ~13 M_⊗, is located at 50 au. In both the low- and moderate-turbulence disks (α_diff = 10⁻⁵ and 10⁻⁴), the Neptunelike planet can generate the dust ring and gap whose widths are a few au in the distribution of ~ 0.2 mm sized dust within 1 Myr.

5.2 Implications for planet formation via pebble accretion

Once the gap forms in the distribution of pebbles, it reduces the accretion rate of pebbles onto protoplanets. As shown in Fig. 17a, even planets with masses below the pebble isolation mass can affect significantly Σ_d(t). When the planetary mass exceeds m ≳ 0.1 (M_p ≳ 0.7 M_⊕ at 1 au), the subsequent growth of the planet would be suppressed and planets remain small in the terrestrial planet forming-region. The suppression of pebble accretion in the terrestrial planet forming-region due to the dust-gap-opening effect by the gas-flow mechanism may be helpful in explaining the current observed period-mass diagram of exo-planets, in which a large fraction of low-mass planets (≲10 M_⊕) has been found at short-period orbits (≲100 days; e.g., Fressin et al. 2013; Weiss & Marcy 2014, see Sect. 4.4.2 of Paper I for further discussion).

5.3 Comparison to disk observations

We compared our semi-analytic models of the widths of the dust ring and gap with the observational data, finding that a fraction of the observed dust rings and gaps could be explained by the gas-flow driven by low-mass planets. We considered here a single planet embedded in a disk. Provided that the thermal emission of the dust is optically thin, and the opacity and the temperature are constant within a dust substructure, the dust surface density is proportional to the observed intensity profile, Σ_d ∝ I_v. In the following paragraphs, we only compared our semi-analytic models with the observed widths of the dust ring and gap. A direct comparison of our semi-analytic model of the dust gap depth with the observational data is difficult because the optically thin assumption would be invalid at the dust ring locations (Guerra-Alvarado et al. 2024; Ribas et al. 2024) and, consequently, the observed intensity does not correspond to a unique dust surface density. In the highly optically thick regime the rings may not be observed, because strong optical depth effects lead to flat-topped radial intensity profiles (Dullemond et al. 2018). Our semi-analytic model of the dust ring width is valid as long as a ring with a moderate optical depth is detectable in the radial intensity profile. The uncertainty in the optical depth of the observed no-flat-topped ring does not significantly affect the observed value of the dust ring width (Dullemond et al. 2018).

To compare with the observational data, we converted the units of ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{SA}}$ and ${\mathcal{W}}_{\text{ring\hspace{0.17em}}}^{\text{SA}}$ from H to au using Eq. (A.3). Same as in Sect. 5.1, we considered the typical steady accretion disk model with the fixed stellar mass, stellar luminosity, and the mass accretion rate (Appendix A): M_* = 1 M_⊙, L_* = 1 L_⊙, and Ṁ_* = 10⁻⁸ M_⊙/yr. We note that rings and gaps have been observed around various types of stars, so that, in reality, these values vary in disks (Huang et al. 2018; Long et al. 2018; Bae et al. 2023): M_* ~ 0.2–2 M_⊙, L_* ~ 0.1–20 L_⊙, and Ṁ_* ~ 10⁻¹⁰–10⁻⁷ M_⊙/yr.

5.3.1 Dust gap width

Figure 19 compares our semi-analytic model of the dust gap width, ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{SA}}$, with the observational data, ${\mathcal{W}}_{\text{gap}}^{\text{obs}}$, in which the dust gap widths are plotted as a function of the gap location, r_gap. The observational data were obtained from the Atacama Large Millimeter/submillimeter Array (ALMA) surveys (Huang et al. 2018; Long et al. 2018; Zhang et al. 2023), including the Disk Substructures at High Angular Resolution Project (DSHARP; Huang et al. 2018). The assumed Stokes number and the Mach number of the headwind to plot ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{SA}}$ were the same values as in the fiducial case: St = 10⁻³, and ℳ_hw = 0.03.

About 20% of the observed dust gaps, whose widths are comparable to the gas scale height ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{obs\hspace{0.17em}}}\sim H$, could be explained by our gas-flow mechanism in the moderate-turbulence disks (α_diff = 10⁻⁴; Fig. 19a). The gas-flow mechanism has the potential to explain the observed wide dust gaps with ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{obs\hspace{0.17em}}}\gtrsim H$ by assuming weaker midplane turbulence (α_diff = 10⁻⁵) and a long time evolution exceeding t ≳ 10⁴. Given that an upper limit for the time required for gap formation is 3 Myr, ~65% of the observed gaps can be explained (Fig. 19b).

These comparisons may suggest the existence of low-mass planets (m ≳ 0.05) at wide orbits as an origin of the observed dust gaps, which could be consistent with a large population of such planets inferred from a population synthesis model (Drazkowska et al. 2023). However, it is difficult to constrain the masses of unseen planets, because the dust gap widths in our model converge when m ≳ 0.3 or t ≳ 10⁴ (Fig. 12).

Figure 19 suggests that low-mass planets within ≲ 10 au could carve gaps as wide as their location, which are entering the transition disk regime (Francis & van der Marel 2020). However, since the dependence of the evolution time of the dust cavities on the parameters (m, ℳ_hw, St, and α_diff) is unclear (Fig. 6), further investigations are needed to link the transition disks with our gas-flow mechanism.

Fig. 19 Dust gap width as a function of the dust gap location. We set St = 10−3 and ℳhw = 0.03, and varied the turbulent parameter in each panel, αdiff = 10−4 (panel a) and αdiff = 10−5 (panel b). The blue and orange solid lines denote H and 5 H, respectively (Eq. (A.3)). Panel a: Blue shaded region given by ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{SA}}\left(t\right)$ (Eq. (26)) at t ≤ 104. The lower and upper limits were set by m = 0.05 and 0.5, respectively. Panel b: Dashed lines given by Eq. (26) with a fixed planetary mass, m = 0.1 . The different colors correspond to different times, t ≤ 104, t = 3 × 104, 105, and 3 × 105, respectively. We hatched the region in which the time required for gap formation exceeds 3 Myr for the assumed dimensionless time, t. The observational data indicated with pink and purple markers show the observed gaps that are accompanied by a ring, taken from Zhang et al. (2023) (compiled from Huang et al. (2018); Long et al. (2018), and Zhang et al. (2023)) and Huang et al. (2018) (DSHARP sample), respectively. These samples do not include disks with inner dust cavities.

Fig. 20 Dust ring width as a function of the ring location. We set ℳhw = 0.03 and αdiff = 10−4. The dashed and dotted lines are given by Eq. (33) with St = 10−3 and 10−4, where we converted the units of ${\mathcal{W}}_{\text{ring\hspace{0.17em}}}^{\text{SA}}$ from H to au using Eq. (A.3). The blue and orange solid lines denote H and 5 H, respectively. We hatched the region where the time required for ring formation exceeds 3 Myr for the assumed dimensionless time, t. The observational data with pink and purple markers are from Bae et al. (2023) and Huang et al. (2018) (DSHARP sample), respectively.

5.3.2 Dust ring width

The observed dust ring widths range from a few au to a few tens of au, which are predominantly wider than those predicted by our semi-analytic model. Figure 20 compares our semi-analytic model of the dust ring width with the observational data, ${\mathcal{W}}_{\text{ring\hspace{0.17em}}}^{\text{obs\hspace{0.17em}}}$, in which the dust ring widths are plotted as a function of the ring locations. The observational data were compiled from Huang et al. (2018) and Bae et al. (2023). The assumed Mach number of the headwind and the turbulent parameter to plot ${\mathcal{W}}_{\text{ring\hspace{0.17em}}}^{\text{SA}}$ were the same values as in the fiducial case: ℳ_hw = 0.03 and α_diff = 10⁻⁴. We considered St = 10⁻⁴ in Fig. 20. We note that we do not show α_diff = 10⁻⁵ or St = 10⁻³ because it leads to narrower gaps and poor agreement with observed rings.

About 15% of the observed rings, whose widths are less than ${\mathcal{W}}_{\text{ring\hspace{0.17em}}}^{\text{obs\hspace{0.17em}}}\lesssim H$, could be explained by the gas-flow mechanism, given that an upper limit of the time required for ring formation is 3 Myr (Fig. 20). We found that only ~3% of the observed dust rings could be explained by the gas-flow mechanism when St = 10⁻³ and α_diff = 10⁻⁴. This may suggest that the dust particles at rings are small due to collisional fragmentation or bouncing (Blum & Wurm 2008; Güttler et al. 2010; Zsom et al. 2010). The inferred Stokes numbers at the dust rings in this study are consistent with the results of dust growth simulations considering the fragility of porous dust (St ~ 10⁻⁴–10⁻³; Ueda et al. 2024). However, when compact dust was considered, larger Stokes numbers were inferred from the multi-wavelength analysis of the continuum emission, the modeling of dust rings at gas pressure maxima, and the dust growth simulations (St ~ 10⁻³–10⁻¹; Sierra et al. 2021; Doi & Kataoka 2023; Jiang et al. 2024). Further discussion is difficult, because the observed wide rings might not sufficiently resolved, and thus they could be composed of multiple narrow rings (Bae et al. 2023).

5.3.3 Potential of multiple planets

So far we have been only considered the dust ring and gap formation by a single planet. The observed wide dust gaps with ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{obs\hspace{0.17em}}}\gtrsim H$ may also be explained by the radially-outward gas flows induced by the multiple planets. When the multiple planets are embedded in the disks with an orbital separation of $\text{Δ}a\lesssim {\mathcal{W}}_{\text{gap}}^{\text{SA}}$, the wide dust gaps may form which are shared by the multiple planets, although the orbital stability of planets is beyond the scope of this study. In this case, a single dust ring forms outside the orbit of the outermost planet.

The observed wide rings may consist of the narrow rings. When the multiple planets are embedded in the disks with the orbital separation of $\text{Δ}a>{\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{SA}}\sim H$, the multiple rings may form with a separation of ~∆a. If the spatial resolution of the observations is low (≳ H), these multiple narrow rings may be observed as a single wide ring.

5.4 Implications for future disk observations

The relatively deep and wide dust gaps formed by the gas flows driven by low-mass planets in the process of their formation in the inner few au of low-turbulence disks could be detected by future observations (Fig. 17a). Due to limitations of the angular resolution, it is difficult to detect the dust substructures in the inner few au of disks by current observations. A possible future extension would improve the angular resolution of the current ALMA by several times, which could lead to further detection of dust substructure in the inner few au of disks (Carpenter et al. 2020; Burrill et al. 2022). A next-generation Very Large Array (ngVLA) is expected to capture the dust thermal emission at ~mm–cm wavelengths with the best angular resolution of ~0.001 arcsec (Selina et al. 2018), which will also provide the capability to probe the inner few au of disks. Simulations of ngVLA observations suggest that the dust gaps with the widths of ~2–3 au are detectable at ~5 au under weak turbulence (α_diff ≲ 10⁻⁵; Ricci et al. 2018; Harter et al. 2020).

Future high angular resolution observations may also detect narrow dust rings and gaps with widths comparable to or less than the gas scale height, H. Thus, our gas-flow mechanism needs to be compared with these future observations.

Future observations may detect the outflow region produced by our gas-flow mechanism. The spatial scale of the outflow region is on the order of the disk gas scale height, ~ H ≈ 8.9 au(r/100au)^9/7 (Eq. (A.4)), in the radial and azimuthal directions (Kuwahara & Kurokawa 2024). The maximum amplitude of the velocity perturbation of the outflow is ~0.3 c_s = 0.03 υ_K(h/0.1)≃0.09km/s(r/100au)^−1/2(M_*/M_⊙)^1/2 (Kuwahara & Kurokawa 2024). Given that the distance to the disk is 100 parsecs and a planet is embedded in the disk at 100 au, the capability to resolve the disk at ~ 0.1 arcsec and ~ 0.1 km/s is required. The current ALMA has the angular resolution of ≳ 0.1 arcsec and the velocity resolution of ≳ 0.01 km/s for the gas. The kinematic features of the outflow induced by low-mass planets may be detectable. However, as discussed in Sect. 4.5.2 of Paper I, the kinematic features of the outflow can only be appeared in the region close to the midplane (ɀ < H; Fig. 3). Therefore, molecules that can trace low heights in disks, such as C¹7O, HCN, and C₂H, should be used as tracers.

5.5 Comparison to previous studies

Previous studies have mostly focused on gap-opening planets to explain the observed dust gap widths (the gas-gap mechanism), often using empirical relations between the planetary mass and the dust gap width obtained from the disk-planet interaction simulations. Zhang et al. (2018) performed 2D hydrodynamical simulations of gas and dust with gas-gap-opening planets. The authors defined the dust gap width by ∆^Z18 ≡ (r_out – r_in)/r_out, where r_in and r_out are the edges of the dust gap normalized by the planet location r_p, and then obtained ∆^Z18 ~ 0.1–1 for different disk models with h = 0.05, 0.07, and 0.1. In our dimensionless unit, the dust gap width defined in Zhang et al. (2018) can be described by ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{Z18}={\text{Δ}}^{\text{Z}18}{r}_{\text{out\hspace{0.17em}}}/h$. Assuming r_out ~ r_p ≡ 1, we obtained ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{Z18}~1-20H$.

An empirical relation in which the dust gap width is assumed to be proportional to the Hill radius has been obtained by the hydrodynamical simulations with gas-gap-opening planets (${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{Hill\hspace{0.17em}}}~4-7.5R_{\text{Hill\hspace{0.17em}}}$; Rosotti et al. 2016; Lodato et al. 2019; Wang et al. 2021). In our parameter sets, the Hill radius ranges from R_Hill ≃ 0.22 to 0.55, which leads to ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{Hill\hspace{0.17em}}}~1-4H$.

Dong et al. (2018) performed 2D hydrodynamical simulations of gas and dust with embedded planets, investigating the dust gap formation due to the shallow gas-gap opening by a planet under the weak turbulence condition (α_diff ≲ 10). The authors considered the planets with m = 0.04–0.8, obtaining the empirical relations between the dust gap width and the planetary mass, ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{D}18}\simeq 3.6l_{\text{sh}}$, where l_sh is the so-called shocking length of the density waves launched by an embedded planet (Goodman & Rafikov 2001): $${l_{\text{sh}}\simeq 0.8{\left(\frac{\gamma +1}{12/5}m\right)}^{-2/5}H\simeq 2{\left(\frac{m}{0.1}\right)}^{-2/5}|}_{\gamma =1.4}H.$$(36)

Here γ = 1.4 is the adiabatic index. Thus, in our dimensionless unit, we obtained ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}^{\text{D}18}~3-10\text{\hspace{0.17em}H}$.

The dust gap widths obtained in this study are narrower than those obtained in the previous studies mentioned above as long as a temporally constant dust gap was considered, ${\mathcal{W}}_{\text{gap\hspace{0.17em}}}\lesssim 4/3H$. The difference in dust gap widths obtained in the previous studies and this study is attributed to different physics and the different parameter range of the Stokes number assumed in each study. The models for dust ring and gap formation driven by higher-mass planets carving gas gaps (the gas-gap mechanisms) are more susceptible to occur when larger Stokes numbers are assumed (St ≳ 10⁻³–10⁻²; Zhu et al. 2012, 2014; Rosotti et al. 2016; Weber et al. 2018). Then the dust particles are trapped at x ≳ H. On the other hand, our gas-flow mechanism works well for smaller solids with St ≲ 10⁻². Such dust particles are trapped by the outflow of the gas at x ≲ 2/3 H.

Although we did not consider the (shallow) gas-gap formation in this study, we would expect that the gas-flow mechanism can coexist with the gas-gap-opening mechanism. The gas-gap mechanism generates the dust gaps with 𝒲_gap ≳ H in the distribution of large dust (St ≳ 10⁻³–10⁻²). Because the small dust particles can pass through the gas gap due to diffusion or the viscous accretion flow (e.g., Rice et al. 2006), the dust gaps also form in the distribution of small dust (St ≲ 10⁻²) by the gas-flow mechanism, whose widths depend on the assumed parameters. The locations of the outer edge of the dust gaps should differ between the gas-flow and gas-gap-opening mechanisms.

Several studies have shown that low-mass planets, which do not form gas gaps or pressure bumps, can create rings and gaps. The gravitational torque exerted by embedded planets can carve gaps that appear only in the dust distribution (Muto & Inutsuka 2009; Dipierro & Laibe 2017). Muto & Inutsuka (2009) showed that, when the global pressure gradient of the disk gas is neglected, a planet with a mass of 2 M_⊕ can open a gap in the distribution of large dust grains (St ≳ 0.1). Dipierro & Laibe (2017) derived a grain-size-dependent criterion for dust gap opening in disks, showing that a planet with a mass of ≲10 M_⊕ can carve a dust gap when St ≳ 1.

Compared to these previous studies that work well for the dust with large Stokes numbers, our results show the opposite trend. In our study, dust rings and gaps form when the smaller Stokes numbers are considered (St ≲ 10⁻²), as we focused on the potential effect of the gas flow driven by an embedded planet, which was not considered in the previous studies.

5.6 Model limitations

Although we did not consider the evolution of the gas disk, the outer part of a disk evolving with viscous angular momentum transport could spread outward due to the conservation of the angular momentum (Lynden-Bell & Pringle 1974). The disk gas at r > r_exp expands outward, where r_exp = r₀/2(1 + t/t_v), r₀ is an initial disk radius, and t_v is a characteristic viscous timescale at r = r₀. The small dust considered in this study could move outward together with the resulting outward flow of the background disk gas (Liu et al. 2022). When a planet is embedded at r > r_exp, the dust particles will be trapped by the radially-inward outflow of the gas induced by an embedded planet. Then the dust ring could form inside the planetary orbit and the dust could deplete outside the planetary orbit. However, the direction of the flow of the background disk gas at the midplane is a controversial issue. The disk gas evolution may be driven by magnetic disk winds (e.g., Bai & Stone 2013). The winds extract angular momentum from the disk surface and drive inward gas accretion at the midplane.

Low-mass planets would undergo inward migration (so-called type I migration; Ward 1986; Tanaka et al. 2002). The timescale of the type I migration is described by $$t_{\text{typeI\hspace{0.17em}}}\simeq \frac{M_{*}}{M_{\text{p}}}\frac{M_{*}}{{\displaystyle {\sum }_{\text{g}}{r}^2}}{h}^3{\text{Ω}}^{-1}={10}^4{\left(\frac{m}{0.1}\right)}^{-1}{\left(\frac{{\text{Σ}}_{\text{g}}{r}^2/M_{*}}{{10}^{-3}}\right)}^{-1}{\text{Ω}}^{-1},$$(37)

where Σ_g is the gas surface density. The type I migration timescale is comparable or shorter than the timescale for the dust ring and gap formation by the gas-flow mechanism. When the migration timescale is shorter than the timescale for the dust ring and gap formation, the positions of the dust ring and gap do not necessarily coincide with the planet location (Kanagawa et al. 2021).

We fixed the planetary masses in this study, whereas planets grow by pebble accretion. A perturbation to the dust surface density becomes strong as the planetary mass increases. The dust gap widens and deepens as the planetary mass increases if the growth timescale of the planet is shorter than the timescale for the dust gap formation. When the planetary mass exceeds approximately m ~ 0.1–0.3, the gap width is independent of the planetary mass (Fig. 12). Thus, we could constrain the lower bound of the mass of the embedded planet from the observed dust gap width. We also fixed the Stokes number, even though it could be varied by dust growth and fragmentation at the dust rings. These additional processes would further complicate the time evolution of Σ_d(t) and should be included in further studies.

An eccentricity of an embedded planet could alter our results. A planet on eccentric orbit induces a time-dependent gas flow field (Bailey et al. 2021), in which the radially-outward gas flow would disappear. However, the eccentricity of the planet can be damped by the disk-planet interaction. The eccentricity damping timescale is given by (Tanaka & Ward 2004) $$\begin{array}{l}{t}_{\text{damp\hspace{0.17em}}}\simeq 1.282\frac{M_{\ast }}{M_{\text{p}}}\frac{M_{\ast }}{{\text{Σ}}_{\text{g}}{r}^2}{h}^4{\text{Ω}}^{-1}\hfill \\ \simeq 6.4\times {10}^2{\left(\frac{m}{0.1}\right)}^{-1}{\left(\frac{{\text{Σ}}_{\text{g}}{r}^2/M_{\ast }}{{10}^{-3}}\right)}^{-1}\left(\frac{h}{0.05}\right){\text{Ω}}^{-1}.\hfill \end{array}$$(38)

Thus, the dust ring and gap formation by the gas-flow mechanism is valid when the eccentricity damping timescale is shorter than the timescale for the dust ring and gap formation.

We finally note that the backreaction of dust on gas, which is not considered in this study, would be important at the dust rings. When the backreaction is included, the axisymmetric dust rings form without planets due to the self-induced dust trap mechanism (Gonzalez et al. 2017; Vericel & Gonzalez 2020; Vericel et al. 2021). When a planet embedded in a disk, the dust ring outside the planetary orbit with the high local dust-to-gas density ratio (≳1) could be unstable and broken into small-scale dust-gas vortices (Pierens et al. 2019; Yang & Zhu 2020), which could change an axisymmetric morphology of the dust rings considered in this study.

6 Conclusions

We investigated the time evolution of dust rings and gaps formed by low-mass planets inducing a radially-outward gas flow. By fitting our numerical results, we developed semi-analytic models describing the widths of the dust rings and gaps and the depth of the dust gaps. Our main results are as follows:

1. Under weak turbulence (α_diff ≲ 10⁻⁴), low-mass planets with m > 0.05 (corresponding to >0.33 M_⊕ at 1 au or ≳1.7 M_⊕ at 10 au) can generate dust rings and gaps in the distribution of small dust, St ≲ 10⁻².

2. Dust gaps have a width comparable to the gas scale height, but can expand further in size when m ≳ 0.1 and α_diff ≲ 10⁻⁵, at a rate set by the dust drift speed (Eq. (26)).

3. The dust gap depth deepens as the planetary mass increases when m ≲ 0.3, but converges at m ≳ 0.3 to a depletion factor of δ_gap ~ 0.2 when α_diff = 10⁻⁴ (δ_gap ~ 10⁻⁷ when α_diff = 10⁻⁵; Eq. (28)). Deeper dust gaps form when smaller turbulent parameters are assumed.

4. The dust rings outside of the planetary orbit widen with time due to diffusion and then reach a steady state. The widths range from ~ 0.1 H to 10 H depending on St, α_diff, and ℳ_hw (Eq. (33)).

By comparing our semi-analytic models of the dust rings and gaps with the observational data, we found that up to approximately 65% (15%) of the observed dust gaps (rings) could be generated by the gas-flow driven by a single low-mass planet. When St = 10⁻³ and α_diff = 10⁻⁴ are considered as the fiducial values, low-mass planets could explain approximately 20% (3%) of the observed dust gaps (rings) with the radial widths of ~ H within t ≤ 10⁴-10⁵, corresponding to ≤ 0.05–0.5 Myr at 10 au. On longer timescales (t ≳ 10⁴–10⁵), the gas-flow mechanism also has the potential to explain approximately 65% (15%) of the observed wide gaps (rings) with widths exceeding the gas scale height H. Wide gaps require a low level of midplane turbulence (α_diff ≲ 10⁻⁵) and wide rings require very small Stokes numbers (St ≲ 10⁻⁴).

Our model for dust ring and gap formation favors low values of St (St ≲ 10⁻⁴–10⁻³), which may suggest the existence of fragile dust grains in protoplanetary disks (Okuzumi & Tazaki 2019; Jiang et al. 2024; Ueda et al. 2024). A fraction of the observed disk substructures may already be consistent with such low-mass planets in wide orbits where they may be ubiquitous during planet formation, before migrating into their final orbits (Drazkowska et al. 2023).

Data availability

The additional figures are available on Zenodo: https://doi.org/10.5281/zenodo.13982345

Acknowledgements

We would like to thank an anonymous referee for constructive comments that have improved the quality of this manuscript. We thank Athena++ developers. Numerical computations were in part carried out on Cray XC50 at the Center for Computational Astrophysics at the National Astronomical Observatory of Japan. A.K. and M.L. acknowledge the ERC starting grant 101041466-EXODOSS. H.K. is supported by JSPS KAKENHI Grant Nos. 20KK0080, 21H04514, 21K13976, 22H01290, 22H05150.

Appendix A Conversion to dimensional quantities

It is practical to convert dimensionless quantities into dimensional ones for discussion. The following sections describe the method of the conversion for a given disk model.

Appendix A.1 Disk model

We considered the typical steady accretion disk model with a dimensionless viscous alpha parameter (Shakura & Sunyaev 1973), α_acc, including viscous heating due to the gas accretion and stellar irradiation heating (e.g., Ida et al. 2016). For simplicity, we fixed the stellar mass, the stellar luminosity, the mass accretion rate, and the viscous alpha parameter as M_* = 1M_⊙, L_* = 1L_⊙, Ṁ_* = 10⁻⁸M_⊙/yr, and α_acc = 10⁻³.

The disk midplane temperature is given by T_disk = max(T_vis, T_irr), where T_vis and T_irr are temperatures determined by viscous heating and stellar irradiation (Garaud & Lin 2007; Oka et al. 2011; Ida et al. 2016), $$T_{\text{vis}}\simeq 200{\left(\frac{M_{*}}{1M_{\odot }}\right)}^{3/10}{\left(\frac{{\alpha }_{\text{acc}}}{{10}^{-3}}\right)}^{-1/5}{\left(\frac{{\dot{M}}_{*}}{{10}^{-8}{M}_{\odot }/\text{yr}}\right)}^{2/5}{\left(\frac{r}{1\text{\hspace{0.17em}au}}\right)}^{-9/10}\text{\hspace{0.17em}K},$$(A.1) $$T_{\text{irr}}\simeq 150{\left(\frac{L_{*}}{1L_{\odot }}\right)}^{2/7}{\left(\frac{M_{*}}{1M_{\odot }}\right)}^{-1/7}{\left(\frac{r}{1\text{\hspace{0.17em}au}}\right)}^{-3/7}\text{\hspace{0.17em}K}.$$(A.2)

The disk gas scale height is given by: $$H=\frac{c_{\text{s}}}{\text{Ω}}=\sqrt{\frac{k_{\text{B}}{T}_{\text{disk}}}{\mu m_{\text{p}}}}\frac{1}{\text{Ω}},$$(A.3)

where k_B is the Boltzmann constant, µ = 2.34 is the mean molecular weight, and m_p is the proton mass. Thus, the aspect ratio of the disk is given by $$h\equiv \frac{H}{r}=\max \left(h_{\text{g},\text{\hspace{0.17em}vis\hspace{0.17em}}},h_{\text{g},\text{irr}}\right),$$(A.4)

where $$h_{\text{g},\text{vis}}\simeq 0.027{\left(\frac{M_{*}}{M_{\odot }}\right)}^{-7/20}{\left(\frac{{\alpha }_{\text{acc}}}{{10}^{-3}}\right)}^{-1/10}{\left(\frac{{\dot{M}}_{*}}{{10}^{-8}{M}_{\odot }/\text{yr}}\right)}^{1/5}{\left(\frac{r}{1\text{\hspace{0.17em}au}}\right)}^{1/20},$$(A.5) $$h_{\text{g},\text{irr}}\simeq 0.024{\left(\frac{L_{*}}{L_{\odot }}\right)}^{1/7}{\left(\frac{M_{*}}{M_{\odot }}\right)}^{-4/7}{\left(\frac{r}{1\text{\hspace{0.17em}au}}\right)}^{2/7}.$$(A.6)

The gas surface density is given by: $${\text{Σ}}_{\text{g}}=\frac{{\dot{M}}_{*}}{3\pi {\alpha }_{\text{acc}}{H}^2{\text{Ω}}_{\text{K}}}=\min \left({\text{Σ}}_{\text{g},\text{vis}},{\text{Σ}}_{\text{g},\text{irr}}\right),$$(A.7)

where $${\text{Σ}}_{\text{g},\text{vis}}\simeq 2.1\times {10}^3\text{\hspace{0.17em}g}/{\text{cm}}^2{\left(\frac{0.027}{h}\right)}^{-2}{\left(\frac{M_{*}}{M_{\odot }}\right)}^{1/5}{\left(\frac{{\alpha }_{\text{acc}}}{{10}^{-3}}\right)}^{-4/5}{\left(\frac{\dot{M}*}{{10}^{-8}M\odot /\text{yr}}\right)}^{3/5}{\left(\frac{r}{1\text{\hspace{0.17em}au}}\right)}^{-3/5},$$(A.8) $${\text{Σ}}_{\text{g},\text{irr}}\simeq 2.7\times {10}^3\text{\hspace{0.17em}g}/{\text{cm}}^2{\left(\frac{L_{*}}{L_{\odot }}\right)}^{-2/7}{\left(\frac{M_{*}}{M_{\odot }}\right)}^{9/14}{\left(\frac{{\alpha }_{\text{acc}}}{{10}^{-3}}\right)}^{-1}\left(\frac{\dot{M}*}{{10}^{-8}M\odot /\text{yr}}\right){\left(\frac{r}{1\text{\hspace{0.17em}au}}\right)}^{-15/14}.$$(A.9)

Fig. A.1 Mach number, planetary mass, and the dust size as a function of the orbital radius. The gray shaded region in panel c denotes the Stokes regime.

Appendix A.2 Orbital radius of the planet

When we convert a dimensionless quantity into a dimensionless one, we need to determine the orbital distance of the planet, r_p. To do this, we specify r_p based on the value of the Mach number of the headwind, ℳ_hw.

From Eq. (A.4), the Mach number of the headwind is given by: $${ℳ}_{\text{hw}}=-\frac{h}{2}\frac{\text{\hspace{0.17em}d\hspace{0.17em}}\ln p}{\text{\hspace{0.17em}d\hspace{0.17em}}\ln r}\simeq \max \left(0.034{\left(\frac{r_{\text{p}}}{1\text{\hspace{0.17em}au}}\right)}^{1/20},0.033{\left(\frac{r_{\text{p}}}{1\text{\hspace{0.17em}au}}\right)}^{2/7}\right),$$(A.10)

where we set d ln p/d ln r = −2.55 for the viscous heating regime and dln p/d ln r = −2.78 for the irradiation heating regime, respectively (Ida et al. 2016). From Eq. (A.10), we can set r_p ≃ 1 au when ℳ_hw = 0.03 (r_p ≃ 50 au when ℳ_hw = 0.1) as a reference value of the planet location (Fig. A.1a).

Appendix A.3 Planet mass

In the disk model given in Sect. A.1, the planet mass is calculated by: $$M_{\text{p}}=m{M}_{*}{h}^3,$$(A.11) $$\simeq \max \left(6.6m{\left(\frac{r_{\text{p}}}{1\text{\hspace{0.17em}au}}\right)}^{3/20},4.6m{\left(\frac{r_{\text{p}}}{1\text{\hspace{0.17em}au}}\right)}^{6/7}\right)M_{\oplus }$$(A.12) $$\simeq \{\begin{array}{l}0.66\left(\frac{m}{0.1}\right)M_{\oplus }\left(r_{\text{p}}=1\text{\hspace{0.17em}au}\right),\hfill \\ 13\left(\frac{m}{0.1}\right)M_{\oplus }\left(r_{\text{p}}=50\text{\hspace{0.17em}au}\right).\hfill \end{array}$$(A.13)

The pebble isolation mass is given by (Bitsch et al. 2018): $$M_{\text{iso}}=25M_{\oplus }{\left(\frac{h}{0.05}\right)}^3\left[0.34{\left(\frac{3}{\log {\alpha }_{\text{diff\hspace{0.17em}}}}\right)}^4+0.66\right].$$(A.14)

Figure A.1b shows M_p for different m as a function of the orbital radius.

Appendix A.4 Dust size

From Eq. (4), the physical size of dust, s, is calculated by: $$s=\min \left(\frac{{\text{Σ}}_{\text{g}}\text{St}}{\sqrt{2\pi }{\rho }_{•}},{\left(\frac{9\mu m_{\text{H}}H\text{St}}{4{\rho }_{•}{\sigma }_{\text{mol}}}\right)}^{1/2}\right).$$(A.15)

In Eq. (A.15), we used the following equations for the stopping time of dust: $$t_{\text{stop\hspace{0.17em}}}=\{\begin{array}{l}\frac{\sqrt{2\pi }{\rho }_{•}s}{{\text{Σ}}_{\text{g}}\text{Ω}}(s\le 9\lambda /4:\text{\hspace{0.17em}Epstein\hspace{0.17em}regime\hspace{0.17em}}),\hfill \\ \frac{4{\rho }_{•}{\sigma }_{\text{mol}}{s}^2}{9\mu m_{\text{p}}H\text{Ω}}(s>9\lambda /4:\text{\hspace{0.17em}Stokes\hspace{0.17em}regime\hspace{0.17em}}),\hfill \end{array}$$(A.16)–(A.17)

where ρ_• = 2 g/cm³ is the internal density of dust and σ_mol = 2 × 10⁻¹⁵ cm² is the molecular collision cross section.

In our disk model, the gas drag regime switches from the Stokes to the Epstein regime when St ≤ 6.6 × 10⁻³ at r_p = 1 au (St ≤ 2.2 × 10³ at 50 au). In the Epstein regime, we have $$s\simeq \{\begin{array}{l}3.7\text{\hspace{0.17em}mm}\left(\frac{\text{St}}{{10}^{-3}}\right)\left(\frac{{\text{Σ}}_{\text{g}}}{2.1\times {10}^3\text{\hspace{0.17em}g}/{\text{cm}}^2}\right){\left(\frac{{\rho }_{•}}{2\text{\hspace{0.17em}g}/{\text{cm}}^3}\right)}^{-1}\left(r_{\text{p}}=1\text{\hspace{0.17em}au}\right),\hfill \\ 0.076\text{\hspace{0.17em}mm}\left(\frac{\text{St}}{{10}^{-3}}\right)\left(\frac{{\text{Σ}}_{\text{g}}}{41\text{\hspace{0.17em}g}/{\text{cm}}^2}\right){\left(\frac{{\rho }_{•}}{2\text{\hspace{0.17em}g}/{\text{cm}}^3}\right)}^{-1}\left(r_{\text{p}}=50\text{\hspace{0.17em}au}\right).\hfill \end{array}$$(A.18)–(A.19)

Figure A.1c shows the dust size for different St as a function of the orbital radius.

Appendix A.5 Time

The orbital period is given by: $$T_{\text{orb\hspace{0.17em}}}\simeq \left(\frac{t}{2\pi }\right){\left(\frac{r_{\text{p}}}{1\text{\hspace{0.17em}au}}\right)}^{3/2}{\left(\frac{M_{*}}{M_{\odot }}\right)}^{-1/2}\text{yr},$$(A.20)

where t is the dimensionless time in our simulations.

Fig. A.2 Time evolution of the minimum dust surface density for different planetary masses. We fixed the Stokes number and the Mach number St = 10−3 and ℳhw = 0.03. We set αdiff = 10−4 in panel a and αdiff = 10−5 in panel b. The square symbols denote the global minimum of Σd,min(t). The horizontal dashed lines mark Σcrit (Eq. 21).

Appendix A.6 Length scale

For each dust surface density simulation, we used a 1D simulation domain with the constant spatial intervals, Δx = 0.01 H (Sect. 2.4). For a given orbital distance of the planet, r_p, we convert the units of the grid cell from H to au based on the following methods. We defined the radial distance of the i-th grid from the central star in au units as r_i, where i = 0 corresponds to the planet location and i = 1, 2, … (i = −1, −2, …) corresponds to the position outside (inside) the planetary orbit. Assuming that the disk aspect ratio is constant, we computed r_i by the following equation: $$r_0=r_{\text{p}},$$(A.21) $$r_i\simeq \{\begin{array}{ll}{r}_{i-1}\left(1+h\left(r_{\text{p}}\right)\times \text{Δ}x\right)\hfill & (i=1,2,\dots ),\hfill \\ r_{i+1}\left(1-h\left(r_{\text{p}}\right)\times \text{Δ}x\right)\hfill & (i=-1,-2,\dots ).\hfill \end{array}$$(A.22)

Fig. A.3 Global minimum of the dust surface density as a function of the planetary mass. We fixed the Stokes number St = 10−3. We set αdiff = 10−4 in panel a and αdiff = 10−5 in panel b. Different symbols correspond to different Mach numbers. The red solid line is the fitting formula for the numerical results of ℳhw = 0.03 (Eq. B.4). The gray thin lines show the uncertainties of Eq. (B.4).

Fig. A.4 Steady-state dust gap depth as a function of the planetary mass. We fixed the Stokes number St = 10−3. We set αdiff = 10−4 in panel a and αdiff = 10−5 in panel b. Different symbols correspond to different Mach numbers. The red solid line is the fitting formula for the numerical results of ℳhw = 0.03 (Eq. B.9). The gray thin lines show the uncertainties of Eq. (B.9).

Fig. A.5 Steady-state dust ring width as a function of the turbulent parameter. The numerical results for different planetary masses were averaged and plotted with the circle symbols. The solid and dashed lines are given by Eq. (B.12). The shaded regions show the uncertainties.

Fig. A.6 Time-dependent dust ring width as a function of time. The numerical results for different planetary masses were averaged and plotted with the square symbols. The solid and dashed lines are given by Eq. (B.13). The shaded regions show the uncertainties.

Appendix B Fitting formulae for ${\text{Σ}}_{\min }^{\text{fit}}$,${\delta }_{\infty }^{\text{fit}}$, and ${\mathcal{W}}_{\text{ring,}\infty }^{\text{fit}}$

We introduce the fitting formulae for the global minimum of the time-dependent dust surface density, ${\text{Σ}}_{\min }^{\text{fit}}$ the steady-state dust gap depth, ${\delta }_{\infty }^{\text{fit}}$, and the steady-state dust ring width, ${\mathcal{W}}_{\text{ring,}\infty }^{\text{fit}}$. For the fitting processes, we used the numerical results of St ≤ 10⁻³.

Figure A.2 shows the time evolution of the minimum dust surface density for different planetary masses, Σ_d,min(t), obtained from our numerical simulations. The minimum dust surface density has the complex dependence on time. We marked the global minimum of the time-dependent dust surface density, $\underset{t>0}{\min }{\text{Σ}}_{\text{d},\min }(t)$, with the square symbol in Fig. A.2. We considered that when $\underset{t>0}{\min }{\text{Σ}}_{\text{d},\min }(t)<{\text{Σ}}_{\text{crit}}$ the dust gap expands with time (Sect. 4.1).

Figure A.3 shows $\underset{t>0}{\min }{\text{Σ}}_{\text{d},\min }(t)$ as a function of the planetary mass. We found that the dependence of min $\underset{t>0}{\min }{\text{Σ}}_{\text{d},\min }(t)$ on the Mach number is weak. Thus, for the fitting process we used the numerical results of ℳ_hw = 0.03 as the representative value.

We assumed that the numerical result can be fitted by the following sigmoid curve: $${\log }_{10}{\displaystyle {\sum }_{\min }^{\text{fit}}=a\times }\mathrm{erf}\left(\left|b\right|m^c\right).$$(B.1)

Using the scipy.optimize.curve_fit library of python, we constrained the fitting coefficients in Eq. (B.1): a, b, and c. We obtained $$\{\begin{array}{l}a=-0.3713\pm 9.512\times {10}^{-3},\hfill \\ b=321.2\pm 455.3,\hfill \\ c=2.8\pm 0.5927,\hfill \end{array}$$(B.2)

when α_diff = 10⁻⁴ and $$\{\begin{array}{l}a=-4.195\pm 9.930\times {10}^{-2},\hfill \\ b=479.5\pm 469.5,\hfill \\ c=2.8\pm 0.4009,\hfill \end{array}$$(B.3)

when α_diff = 10⁻⁵. From Eqs. (B.2) and (B.3), we obtained $${\log }_{10}{\text{Σ}}_{\min }^{\text{fit}}=-0.37{\left(\frac{{\alpha }_{\text{diff}}}{{10}^{-4}}\right)}^{-1.1}\times \text{erf}\left(3.2\times {10}^2{\left(\frac{{\alpha }_{\text{diff}}}{{10}^{-4}}\right)}^{-0.17}{m}^{2.8}\right),$$(B.4)

where we used the best fit values of the fitting coefficients. Equation (B.1) gives ${\text{Σ}}_{\min }^{\text{fit}}>1$ when α_diff ≳ 10⁻⁴, which is unphysical.

Thus, we set the upper limit and obtained the semi-analytic formula for the global minimum of the dust surface density: $${\text{Σ}}_{\min }^{\text{global}}=\min \left(1,{\text{Σ}}_{\text{min}}^{\text{fit}}\right).$$(B.5)

We plotted Eq. (B.5) in Fig. A.3 with the solid line.

Figure A.4 shows the steady-state dust gap depth as a function of the planetary mass. Similar to the fitting process of we used the numerical results of ℳ_hw = 0.03 as the representative value for the fitting and assumed that the numerical result can be fitted by: $${\log }_{10}{\delta }_{\infty }^{\text{fit}}=a^{\prime }\times \text{erf}\left(\left|b^{\prime }\right|m^{c^{\prime }}\right).$$(B.6)

Using the scipy.optimize.curve_fit library, we obtained the fitting coefficients a′, b′, and c′ as follows: $$\{\begin{array}{l}{a}^{\prime }=-0.6263\pm 2.123\times {10}^{-3},\hfill \\ b^{\prime }=420.6\pm 639.6,\hfill \\ c^{\prime }=2.8\pm 0.6286,\hfill \end{array}$$(B.7)

when α_diff = 10⁻⁴ and $$\{\begin{array}{l}{a}^{\prime }=-7.098\pm 0.2133,\hfill \\ b^{\prime }=399.9\pm 558.9,\hfill \\ c^{\prime }=2.8\pm 0.5792,\hfill \end{array}$$(B.8)

when α_diff = 10⁻⁵. From Eqs. (B.7) and (B.11), we obtained $${\log }_{10}{\delta }_{\infty }^{\text{fit}}=-0.63{\left(\frac{{\alpha }_{\text{diff\hspace{0.17em}}}}{{10}^{-4}}\right)}^{-1.1}\times \mathrm{erf}\left(4.2\times {10}^2{\left(\frac{{\alpha }_{\text{diff\hspace{0.17em}}}}{{10}^{-4}}\right)}^{0.022}{m}^{2.8}\right).$$(B.9)

Equation (B.6) gives ${\delta }_{\infty }^{\text{fit}}>1$ when α_diff ≳ 10⁻⁴. To avoid unphysical solutions, we set the upper limit: $${\delta }_{\infty }=\min \left({\delta }_0,{\delta }_{\infty }^{\text{fit}}\right).$$(B.10)

We plotted Eq. (B.10) in Fig. A.4 with the solid line.

Figure A.5 shows the steady-state dust ring width as a function of α_diff. As shown in Fig. 15, the dust ring width is weakly dependent on the planetary mass. Thus, we averaged the numerical results for different planetary masses and plotted them in Fig. A.5 with the circle symbols. We assumed that the steady-state dust ring width is proportional to the characteristic length in which the drift timescale coincides with the diffusion timescale, ℒ_eq (Eq. 9): ${\mathcal{W}}_{\text{ring,\hspace{0.17em}},\infty }^{\text{fit}}=C\times {ℒ}_{\text{eq}}$. With the least-squares method, we derived the coefficient: $$C=\{\begin{array}{ll}0.63\pm 0.35\hfill & \left({\alpha }_{\text{diff\hspace{0.17em}}}={10}^{-4}\right),\hfill \\ 2.8\pm 0.54\hfill & \left({\alpha }_{\text{diff\hspace{0.17em}}}={10}^{-5}\right),\hfill \end{array}$$(B.11)

and then obtained $${\mathcal{W}}_{\text{ring,\hspace{0.17em}},\infty }^{\text{fit}}=0.63{\left(\frac{{\alpha }_{\text{diff}}}{{10}^{-4}}\right)}^{-0.65}\times {ℒ}_{\text{eq}}.$$(B.12)

We plotted Eq. (B.12) in Fig. A.5 with the solid lines.

Figure A.6 shows the dust ring width as a function of time. Since the dust ring width is weakly dependent on the planetary mass, same as in Fig. A.5, we averaged the numerical results for different planetary masses at each time and then plotted them in Fig. A.6 with the squares symbols. We assumed that the time-dependent dust ring width can be fitted by the following sigmoid curve: $${\mathcal{W}}_{\text{ring\hspace{0.17em}}}^{\text{SA}}(t)={\mathcal{W}}_{\text{ring\hspace{0.17em}},\infty }^{\text{fit}}\left(1-\frac{1}{1+{\left(t/{\tau }_{\text{ring\hspace{0.17em}}}\right)}^q}\right).$$(B.13)

With the least-squares method, we determined the power of (t/τ_ring)^q: q = 0.42 ± 0.045. We plotted Eq. (B.13) in Fig. A.6 with the solid lines.

References