Multi-wavelength Constraints on Dust Dynamics and Size Evolution in Protoplanetary Disk Rings. I. Method

In PPDs, radial pressure gradients induce sub- and super-Keplerian gas rotation, causing headwind and tailwind effects that drive radial dust drift. A crucial quantity for characterizing this effect is the dimensionless stopping time, or Stokes number, ${\rm St}$. In the Epstein regime, the Stokes number is defined as:

where $a$ is the particle’s equivalent radius, $\Sigma_{\rm g}$ is the gas surface density, and $\rho_{\rm s}$ is the bulk density of the dust particle (Weidenschilling, 1977; Birnstiel et al., 2010). The dust size corresponding to ${\rm St}=1$ is given by:

The dust radial drift velocity due to the pressure gradient is given by (Weidenschilling, 1977; Nakagawa et al., 1986):

where $v_{\rm K}$, $c_{\rm s}$, and $P$ represent the Keplerian velocity, sound speed, and gas pressure, respectively. Particles drift towards regions of higher gas pressure, leading to concentration at local pressure maxima. Generally, larger grains exhibit faster radial drift speeds when ${\rm St}\lesssim 1$.

A further important effect is the diffusion of dust grains due to turbulent gas motion. The radial dust particle diffusion coefficient, derived from the Langevin equation for particle motion in a homogeneous isotropic Kolmogorov turbulence spectrum, is (Youdin & Lithwick, 2007):

where $D_{\rm g}=\nu/{\rm Sc}=\alpha c_{\rm s}^{2}/\Omega_{\rm K}{\rm Sc}$ represents the gas diffusion coefficient, $\Omega_{\rm K}$ is the local Keplerian angular velocity, and $\alpha$ is the Shakura-Sunyaev turbulence parameter (Shakura & Sunyaev, 1973). Here, $D_{\rm g}=\nu$ when the Schmidt number ${\rm Sc}=1$. Dust diffusion is more efficient for smaller particles, characterized by smaller ${\rm St}$, which leads to the spreading of dust grains concentrated at the pressure maximum.

Collisions between dust grains lead to size changes. High-energy collisions result in fragmentation (Ormel et al., 2009; Güttler et al., 2010), occurring when the relative velocity of grains during collisions exceeds a critical fragmentation velocity $v_{\rm frag}$. In typical PPD conditions, turbulence provides the dominant source of relative motion for dust grains, with larger grains experiencing higher collision velocities when ${\rm St}\lesssim 1$.

The maximum Stokes number for turbulence-driven fragmentation is given by (Birnstiel et al., 2012):

which corresponds to a critical particle size $a_{\rm frag}$ in the Epstein regime: ¹¹1Note that we have omitted the correction factor $f_{\rm frag}=0.37$ for the fragmentation limit as suggested by Birnstiel et al. (2012).

Dust evolution in PPDs proceeds on several characteristic timescales, each governing a specific physical process.

The first is the grain growth timescale, which characterizes the rate at which dust grains grow. The typical growth timescale is given by (Brauer et al., 2008; Laune et al., 2020; Birnstiel et al., 2012):

where

is the total dust-to-gas surface density ratio. Once the system reaches coagulation-fragmentation equilibrium, the growth timescale reflects the characteristic time for restoring perturbations.

Second, dust grains are continuously transported toward the local pressure maximum or the central star via radial drift. For dust drifting in a Gaussian pressure bump, Eq. (20), the drift timescale simplifies to:

where $w_{\rm p}$ represents the width of the pressure bump, as defined in Eq. (20).

By comparing the growth and drift timescales, the drift limit-the maximum grain size that can be replenished by growth rather than being lost due to radial drift-is determined. For dust grains within pressure bumps, this limit is ²²2For disk particles drifting toward the central star, the corresponding drift limit is given by(Birnstiel et al., 2012):$a_{\rm drift,*}=\frac{2\Sigma_{\rm d}}{\pi\rho_{\rm s}}\frac{v_{\rm K}^{2}}{c_{\rm s}^{2}}\left|\frac{d\ln P}{d\ln r}\right|^{-1}.$:

Expressed in terms of grain size, it follows as:

The final timescale considered is the dust diffusion timescale, which quantifies the time required for dust to spread throughout the disk. In practice, the diffusion timescale depends on grain size. For simplicity, we focus on well-coupled small grains when characterizing it. Since our analysis centers on the dynamical evolution of dust particles near a pressure bump, the diffusion timescale is given by:

To provide a precise quantitative description of various quantities within a pressure bump and to avoid potential confusion, we introduce a standardized notation system before proceeding with the discussion.

Within the 1D dust ring model, the radial profiles of physical quantities in the central region of a pressure bump are approximated as Gaussian functions. The standard deviation of these Gaussian profiles is defined as the corresponding width. The subscript of each width notation corresponds to the respective quantity.

Table 1 provides a summary of the various ring width definitions used in this study. To characterize the dust distribution in both the radial and size domains, by distributing the total mass into logarithmic size bins, the dust size distribution is defined as:

where $\Sigma_{\rm d}(r)$ represents the vertically integrated total dust surface density. The subscript "$\sigma_{\rm d},a$" of $w_{\sigma_{\rm d},a}$ is used to specify the width of the radial profile $\sigma_{\rm d}(a,r)|_{a=a}$ for a given dust population. The same principle also applies to single species dust-to-gas ratio $\epsilon(a,r)$.

The quantity $w\mathrm{{}_{p}}$ is also used as an approximation for the width of the gas surface density profile, $\Sigma_{\rm{g}}$, as adopted in this study and in Dullemond et al. (2018). This approximation holds if the temperature variation across the bump is negligible.

By balancing dust drift and diffusion, while neglecting dust growth and fragmentation, the intrinsic width of dust rings in surface density for a single dust species, $w_{\sigma_{\rm d},a}$, is given by (Dullemond et al., 2018; Morbidelli, 2020).

where $\psi=\sqrt{\alpha/(\mathrm{St_{0}}(a)\cdot\rm{Sc})}$, $\mathrm{St_{0}}$ (which scales with the dust size $a$) is the Stokes number at the bump center, and $\mathrm{Sc}$ is the Schmidt number, typically assumed to be $\mathrm{Sc}=1$.For large particles where $\psi\ll 1$, it follows that $w_{\sigma_{\rm d},a}\simeq w_{\rm p}\psi$, suggesting dust trapping. Conversely, for small particles where $\psi\gg 1$, $w_{\sigma_{\rm d},a}\simeq w_{\rm p}$, indicating that small dust grains are well-coupled with the gas.

It is essential to distinguish between $w_{\epsilon,a}$ and $w_{\sigma_{\rm d},a}$, as both describe the dust distribution from different perspectives. The dust surface density $\sigma_{\mathrm{d}}$ directly characterizes the spatial distribution, whereas the dust-to-gas ratio $\epsilon$ is dynamically more relevant, governing the net dust diffusion flux. These two representations are related by $\sigma_{\rm{d}}({a},r)=\epsilon({a},r)\cdot\Sigma_{\rm{g}}(r)$, leading to the following relation:

For instance, rewriting the drift-diffusion equilibrium solution (Eq. (15)) in terms of the dust-to-gas ratio gives:

In the limit where $w_{\epsilon,a}\to\infty$, the corresponding dust surface density width $w_{\sigma_{\rm d},a}$ asymptotically approaches $w_{\rm p}$, indicating strong dust-gas coupling. Conversely, by combining Eqs. (16) and (17), we can obtain Eq. (15).

The total pressure profile is prescribed as a pressure bump superimposed on a minimized pressure background, expressed as:

where $P_{\rm bump}(r)$ and $P_{\rm disk}(r)$ represent the pressure profiles of the Gaussian bump and the smooth background, respectively. The factor $F(r)$, detailed below, further modifies the background pressure. These components are illustrated in Figure 2.

To facilitate analytical treatment.The function

is introduced to preserve the Gaussian form of the $P_{\rm bump}(r)$ and ensure a well-defined pressure bump width $w_{\rm p}$. We set $f=9.3$ to maintain the pressure bump nearly unchanged.

The pressure bump is initialized with a Gaussian profile:

where $w_{\rm p}$ denotes the width of the pressure bump. The peak mid-plane pressure is given by $P\mathrm{{}_{bump,0}}=\frac{c_{\rm s,0}^{2}\Sigma_{\rm g,0}}{\sqrt{2\pi}H_{\rm p,0}}$, where $H_{\rm p}$ is the pressure scale height, and the subscript "0" indicates values measured at the pressure maximum $r_{0}$.

The smooth background disk pressure is defined as:

where the surface density profile $\Sigma_{\rm disk}(r)$ follows a cutoff power-law profile (Lynden-Bell & Pringle, 1974):

where $A_{\rm p}$ quantifies the contrast between the background disk and the pressure bump, satisfying $P_{\rm bump}(r_{0})=A_{\rm p}P_{\rm disk}(r_{0})$. The parameter $r_{\rm c}$ represents the critical cutoff radius of the disk. The temperature follows $T=T_{0}(r/r_{0})^{-1/2}$, resulting in a radial sound speed profile $c_{\rm s}\propto r^{-1/4}$ (Chiang & Goldreich, 1997).

For consistency, we select fiducial model parameters appropriate for the ring at 100 au in the HD 163296 system, as summarized in Table 2. The total gas pressure profile $P_{\rm total}$ is shown as the dotted-dashed line in Figure 2, demonstrating that the Gaussian bump remains largely intact without significant alteration from the background pressure.

For simplicity, we neglect viscous accretion in the disk and assume the gas pressure profile remains constant over time. In this context, the back-reaction of dust on the gas is not considered in this study.

The dust component is initialized as micron-sized grains, uniformly distributed in the range of $0.5$ to $1\,\mathrm{\mu m}$, with a global dust-to-gas ratio set at $\epsilon\mathrm{{}_{ini}}\equiv\Sigma\mathrm{{}_{d,ini}}/\Sigma_{\rm g}=0.01$ (shown as the dashed line in Figure 2).

Dust grains moving at different velocities undergo collisions, leading to various outcomes. The relative velocities between dust grains arise from turbulence, Brownian motion, vertical settling, and radial/azimuthal drift, with turbulence typically being the dominant contributor. The collision outcomes considered in our simulations include coagulation, fragmentation, and erosion (We adopt the default coagulation model implemented in DustPy, as described in Stammler & Birnstiel (2022) for further details). Fragmentation becomes increasingly likely as the collision velocity exceeds a critical threshold, i.e., the fragmentation velocity $v_{\rm frag}$; otherwise, grains are more likely to adhere, leading to coagulation.

Upon fragmentation, the resulting fragment mass distribution, $n_{\rm frag}(m)\propto m^{-\xi}$, follows a power law with an index of $\xi=11/6$ (Brauer et al., 2008), corresponding to an MRN-like size distribution (Mathis et al., 1977).

The simulation employs a logarithmic mass grid with a mass resolution of seven bins per decade, a configuration extensively tested in previous works for collision modeling (Stammler & Birnstiel, 2022; Ohtsuki et al., 1990; Drążkowska et al., 2014).

The other parameters for the fiducial simulation are based on the 100-au dust ring in HD 163296, as summarized in Table 2. Given a stellar mass of $M_{*}=1.9M_{\odot}$ (Fairlamb et al., 2015) and a pressure bump width of $w_{\mathrm{p}}=23$ au (Rosotti et al., 2020), the temperature at the pressure maximum is set to $T_{0}=12.3~{}\mathrm{K}$ (Guidi et al., 2022). The disk cutoff radius is chosen as $r_{\rm c}=200$ au (Stammler et al., 2019).

The radial domain of our simulations extends from 10 au to 400 au with a grid resolution of $\sim 0.26\ \mathrm{au}$ (corresponding to $0.05H\rm_{p}$) at the bump center. A constant gradient is applied to the inner boundary and a floor value is used for the outer boundary condition. The large radial domain of the disk with an exponential disk profile ensures that the dust dynamics around the ring is not influenced by the boundary condition.

Figure 2 presents the steady-state(5 Myr) total dust-to-gas ratio $\epsilon_{\mathrm{tot}}$ as a function of radius in the disk (the solid line). The total dust-to-gas ratio profile can be described by a Gaussian function

from which we define the intrinsic dust ring width $w_{\epsilon}$ (to distinguish from the width of observed intensity rings $w_{\rm i}$). We can clearly identify significant dust trapping with $w_{\epsilon}<w_{\rm p}$ in the fiducial run, which is also evident from the red dashed line in Figure 4.

Figure 3 shows the steady-state(5 Myr) dust size distribution $\sigma_{\rm d}$. The bump region does not evolve significantly after about $\mathrm{1\ Myr}$. As is shown, dust particles encounter two barriers during their dynamical evolution, one is the drift barrier $a_{\rm drift}$, another is the fragmentation barrier $a_{\mathrm{frag}}$. It is evident that the maximum size at the center of the pressure bump is primarily restricted by $a_{\mathrm{frag}}$, whereas the edges of the bump are constrained by $a_{\rm drift}$.

First, we examine the spatial distribution of the dust ring. Figure 4 shows the measured widths $w_{\epsilon,a}$ for each dust species as a function of $\mathrm{St}_{0}$, where $\mathrm{St}_{0}\equiv\pi a\rho_{\rm s}/2\Sigma_{\rm g,0}$ is the Stokes number at the pressure maximum, which can be interpreted as the grain size. $w_{\epsilon,a}$ is derived by fitting the radial profile of $\epsilon(a,r)$ for each species with a Gaussian profile as in Eq. (23). The dust-to-gas ratio for each dust species is defined as:

For dust rings regulated by coagulation and fragmentation equilibrium, the ring widths $w_{\epsilon,a}$ are nearly identical across dust species with $\mathrm{St_{0}\lesssim 0.37St_{frag}}(r_{0})$ in our simulations, as indicated by the red solid line.

This contrasts with the drift-diffusion scenario, as illustrated by the blue lines and discussed by Dullemond et al. (2018). If dust collisions are neglected (i.e., only drift and diffusion are considered), the width $w_{\epsilon,a}$ becomes size-dependent. Larger, drift-dominated grains are more tightly concentrated, whereas smaller, diffusion-dominated particles are less effectively trapped. Consequently, for well-coupled dust, $w_{\sigma_{\rm d},a}/w_{\rm p}\simeq 1$, which corresponds to $w_{\epsilon,a}/w_{\rm p}\propto\sqrt{1/{\rm St}}>1$, as explained in Section 2.3.

However, grains near the fragmentation limit, $\mathrm{St}_{\rm frag}$, exhibit a narrower width than $w_{\epsilon}$, as shown in Figure 4. This narrowing arises primarily from two factors. First, for even larger pebbles with shorter drift timescales, drift-diffusion becomes the dominant process, and their distribution follows the drift-diffusion equilibrium solution (blue lines in Figure 4). This explains why the red line converges with the blue line for very large pebbles. Second, the fragmentation limit varies with radius, restricting grains near this threshold to a confined region. Specifically, in the central part of the dust ring, the fragmentation limit decreases outward from the ring center (as indicated by the dotted line in Figure 3), allowing $\sim 1\,\mathrm{cm}$ pebbles to exist only near the center. Consequently, these pebbles are more spatially constrained compared to smaller grains. Together, these mechanisms lead to a narrower radial distribution for large pebbles.

Although this effect has a negligible impact on dust dynamics-since the fraction of dust deviating from the common radial distribution is small-it could significantly affect thermal radiation, providing an observable signature that can be used to test the theory, as discussed in the following paper.

Next, we describe the dust size distribution within the central region of the bump. Figure 5 presents the normalized dust size distribution $\sigma_{\rm d}(r)/\Sigma_{\rm d}$ as a function of dust size at different radii represented by the solid (the bump center at 100 and 95 au) and semi-transparent (the wing at 80 and 75 au) lines, respectively. The size distributions at the bump center closely matches with each other, indicating a similar dust size distribution at bump center region, which is notably different from the size distribution in the wings.

The reason for this similar dust size distribution is as follows. The normalized size distribution $\sigma_{\rm d}/\Sigma_{\rm d}$ is primarily characterized by its slope and maximum size. In the case that coagulation-fragmentation is efficient enough that the slope for the dust size distribution is governed by the spatially-independent fragmentation slope $\xi$ (Stammler & Birnstiel, 2022), while the maximum size, determined by turbulence-driven fragmentation, is associated with a spatially invariant Stokes number $\mathrm{St}_{\rm frag}$ (see Eq. 6). Together, these factors suggest that the distribution $\sigma_{\rm d}/\Sigma_{\rm d}$ remains spatially invariant if the dust dynamics are primarily governed by coagulation-fragmentation equilibrium. In contrast, the inclusion of dust drift and diffusion, both of which depend on spatial position, would result in significant variations in $\sigma_{\rm d}/\Sigma_{\rm d}$ within the bump. The similar dust size distribution within the ring thus indicates a coagulation-fragmentation equilibrium throughout the bump center, and suggests that radius-dependent drift and diffusion effects can be neglected in terms of the dust size distribution.

The vertical gray line in Figure 5 represents the predicted fragmentation barrier at the bump center $a_{\rm frag}$, corrected by a factor of $0.37$. The simulation results show that dust accumulates around a slightly smaller size than $a_{\rm frag}$. This discrepancy arises due to the implementation of the Maxwell-Boltzmann distribution for dust-relative velocities (see Eq. 58 in Stammler & Birnstiel, 2022), which is not considered in the analytical analysis. The high-velocity tail of the Maxwell-Boltzmann distribution results in fragmentation occurring at a smaller size.

The size-independent ring width $w_{\epsilon,a}$ aligns with the spatially-insensitive size distribution $\sigma_{\rm d}/\Sigma_{\rm d}$ across all radii, both of which indicate a coagulation-fragmentation equilibrium near the bump center, where the majority of the dust mass resides.

In summary, different dust species within the ring exhibit a similar size and spatial distribution. This suggests that the size distribution at various radii within the ring can be described by the same function, i.e.,

Furthermore, the radial width of the ring for different dust species within the fragmentation limit is represented by a common constant, i.e.,

We skip the discussion about dust size distribution in coagulation-fragmentation equilibrium, $\sigma_{\rm d}$, which is well-studied by Birnstiel et al. (2011).

The width of the dust-to-gas ratio profile for a dust ring, $w_{\epsilon}$ is given by:

The factor of 1.2 is a constant calibrated from simulation results. Here, $\bar{V}$ represents the average drift effect and $\bar{D}$ represents the average diffusion effect for the size distribution at the bump center $\sigma_{\rm d}|_{r_{0}}$, and are computed as

and

A detailed derivation of $w_{\epsilon}$ is provided in Appendix A.

As an illustrative example, consider a power-law size distribution of dust grains, $n(a)\propto a^{-q}$, with a maximum grain size $a\mathrm{{}_{max}}$. If the Stokes number of the largest grain, $\mathrm{St_{max}}$, is less than unity, $\mathrm{St_{max}}<1$, $a_{\mathrm{min}}\ll a_{\mathrm{max}}$ and $q<4$, the following approximation holds:

The ratio $w_{\epsilon}/w_{\rm{p}}$, along with widths for smaller grains of arbitrary size, is determined by $\mathrm{St_{max}}$. In the fragmentation-limited regime, where $\mathrm{St_{max}}=\mathrm{St_{frag}}$, $w_{\epsilon,a<0.37a_{\rm frag}}$ becomes independent of the gas surface density $\Sigma_{\rm g}$, and the dust bulk density $\rho_{\rm s}$ (see Eq. 6). This behavior contrasts with the drift-diffusion scenario, in which the ring width $w_{\epsilon,a<0.37a_{\rm frag}}$ for a given dust species scales as $\sqrt{1/\mathrm{St}(a,\Sigma_{\rm{g}},\rho_{\rm{s}})}$ (Dullemond et al., 2018).

Eq. (25) is verified through comparison with simulation results from an extensive parameter survey. We performed simulations by varying one parameter at a time from the fiducial model and measuring the steady-state intrinsic dust ring width, as given by Eq. (23). The parameters explored include the disk turbulence parameter $\alpha$, the initial dust-to-gas ratio $\epsilon_{\rm ini}$, the initial monomer dust size $a_{0}$, disk temperature $T$, the power-law index of fragments $\xi$, and the fragmentation velocity $v_{\rm frag}$.

We then calculated the ring width using Eq. (25), applying the dust size distributions, $\sigma_{\rm d}|_{r_{0}}$ directly taken from the simulations.

The measured and calculated ring widths are shown in Figure 6. The shaded regions indicate the parameter space where the analytical theory is no longer applicable. Specifically, for very small values of $\alpha$, which correspond to a laminar disk environment, the relative velocity between dust grains never reaches $v_{\rm frag}$, preventing fragmentation. This scenario is out the scope of the current coagulation model (Birnstiel et al., 2009). Additionally, when the dust-to-gas ratio is extremely low, drift and diffusion effects dominate, as discussed in Appendix A (Eqs. A10 and A8), requiring alternative treatment (Appendix A.2).

The ring width calculated using the simulation-derived $\sigma_{\rm d}$ (represented by the dots) matches the measured values perfectly, confirming the validity of the ring width formula described by Eq. (25).

We assume that the dust ring is in drift/diffusion/coagulation/fragmentation equilibrium. And the dust-to-gas ratio is high enough to meet the criterion Eq. (A8). The main steps to reconstruct the dust ring are summarized below, starting from a gas pressure bump:

Step 1

Construct the normalized dust size distribution $(\sigma_{\rm d}/\Sigma_{\rm d})|_{r_{0}}$ at the pressure maximum based on the fragmentation limit $a_{\rm frag}$ and the power-law indices for multi-species dust (see Birnstiel et al., 2018, 2011, for details).

Step 2

Calculate $w_{\epsilon,a}$ for each dust species. First, use the gas bump width $w_{\rm p}$ and the dust size distribution $(\sigma_{\rm d}/\Sigma_{\rm d})|_{r_{0}}$ to determine the dust ring width $w_{\epsilon}$ for the total dust-to-gas ratio $\epsilon$, based on Eq. (25). This width applies to dust species with sizes of ${a<0.37a_{\rm frag}}$. For larger grains, the drift-diffusion balance solution Eq. (17), scaled by the factor $(w_{\epsilon}/w_{\rm p})(\sqrt{\mathrm{St}(0.37a_{\rm frag})/\alpha})$, provides a good approximation.

Step 3

Map the 2D size-spatial distribution of dust, similar to that shown in Figure 3, by distributing $\sigma_{\rm d}(a,r)=\Sigma_{\rm g}(r)\cdot\epsilon(a,r)$ into mass grid bins at various radii.

Appendix A.2 provides an empirical fits for cases when the dust-to-gas ratio is too low to satisfy the criterion Eq. (A8). Utilizing the 2D size-spatial distribution for all dust species, we can reconstruct the dust ring for observational modeling through radiative transfer calculations across various (sub-)mm bands, which will be detailed in our upcoming paper.