© The Authors 2026

1 Introduction

A debris disk can be simply described as a collection of large, approximately kilometer-sized, bodies and the debris produced by their collisions (Krivov & Wyatt 2021). The parent bodies must have formed earlier, in the protoplanetary disk (e.g., Lau et al. 2024), and their mutual collisions result in a steep size distribution that can span at least ten orders of magnitude. The small end of the distribution is set by radiation pressure, which quickly removes the smallest particles from the system (Krivov 2010). At first glance, describing a debris disk requires just a handful of parameters: a birth ring where the parent bodies are located, a size distribution describing the collisional cascade, and the inventory of forces consisting mostly of gravity and radiation pressure (though more refined models can include stellar multiplicity, the presence of planets, Poynting-Robertson drag, self-gravity of the disk, and/or gas drag).

Radiation pressure is a size-dependent force that naturally results in a size segregation as a function of the stellocentric distance (r). The small dust grains are set on highly eccentric orbits (or can become gravitationally unbound to the central star), forming a halo outside the birth ring, while millimeter-sized particles are mostly unaffected and remain close to where they were produced. There is therefore great value in comparing observations taken at different wavelengths. Submillimeter ALMA observations trace the large grains and thus the birth ring, while scattered light observations probe the extended halo as well as the birth ring (Thébault et al. 2023). The ALMA survey to Resolve exoKuiper belt Substructures (ARKS) Large Program (Marino et al. 2026) has provided new submillimeter observations at high angular resolution, informing about the presence of cold CO gas (Mac Manamon et al. 2026) and enabling a detailed comparison with archival Spectro-Polarimetic High contrast imager for Exoplanets REsearch (SPHERE, Beuzit et al. 2019) images. Milli et al. (2026) present a comparison between radial profiles extracted from the ARKS continuum submillimeter observations (see Han et al. 2026) and near-IR scattered light observations. They find that there are six disks with a statistically significant offset, and all but one are also known to host cold CO gas. Conversely, the systems without significant offsets are not known to host gas, except one. As discussed in Takeuchi & Artymowicz (2001), Krivov et al. (2009), Olofsson et al. (2022), and more recently Jankovic et al. (2026) and Weber et al. (2026), gas can alter the dynamics of the dust particles. In this work we therefore investigated whether the presence of gas can explain the presence and strength of the observed offsets. We did not attempt to fit or reproduce any particular set of ALMA and SPHERE observations but rather carried out a parameter study to check under which circumstances such radial offsets can arise.

In Section 2, we describe the simulation setup (from determining the dust density distribution to computing synthetic images). We then present a fiducial model in Section 3 and explore the parameter space to quantify which parameters can explain the radial offsets seen in SPHERE and ALMA observations in Section 4. We discuss our findings in Section 5 before concluding in Section 6.

2 Simulation setup

The end goal of each simulation was to compute the dust density distribution as a function of grain size, accounting for the effects of stellar gravity, radiation pressure, gas drag, and collisions. The approach follows the one originally presented in Thébault (2012), the main differences being the inclusion of gas drag (see also Olofsson et al. 2022) and the fact that unbound grains are included here (similarly to Thébault et al. 2023). In this section we first describe the physical parameters defining our model, then explain how the dust density distribution is determined and how synthetic images are computed.

2.1 Dust and gas parameters

2.1.1 Dust parameters

The full disk can be described using only a few parameters. First, we parametrized the radial distribution of parent bodies in the birth ring by a normal distribution of semimajor axes, centered at a_d with a standard deviation σ_d. Next, the vertical structure was determined by a single parameter ψ, the standard deviation of a normal distribution from which we drew the inclinations of the parent bodies with respect to the midplane.

At the beginning of the simulation, n_part parent bodies are created, which only experience stellar gravity. During the initialization, we drew their mean anomalies from a uniform distribution between 0 and 2π and computed their initial positions and velocities. Each of these “seed” parent bodies release one dust particle (hence n_part particles in total), whose dynamics we followed. Furthermore, a radius (s) was also assigned between a minimum and maximum size, s_min and s_max, respectively, following a power-law distribution dn(s) ∝ s^qds, with a slope q. Finally, each particle was assigned a value of β, the ratio between radiation pressure and gravity, computed following Burns et al. (1979); $$\beta (s)=\frac{3L_{\star }}{16\pi Gc{M}_{\star }}\frac{Q_{\mathrm{p}\mathrm{r}}(s)}{\rho s},$$(1)

where L_* and M_* are the stellar luminosity and mass, G the gravitational constant, c the speed of light, ρ the material density of the dust particles, and Q_pr is the radiation pressure efficiency, computed as $$Q_{\mathrm{p}\mathrm{r}}(s)=Q_{\mathrm{a}\mathrm{b}\mathrm{s}}(s)+Q_{\mathrm{s}\mathrm{c}\mathrm{a}}(s)\times (1-g_{\mathrm{s}\mathrm{c}\mathrm{a}}(s)),$$(2)

averaged over the stellar photospheric model, where Q_abs and Q_sca are the absorption and scattering efficiencies while g_sca is the asymmetry parameter. To have finer control over the blowout size (s_blow below which grains are unbound to the star), we include another free parameter, β_scale, by which we multiply the β(s) curve. For β_scale < 1, this effectively decreases the blow-out size (shifting down β(s)), and vice-versa, for β_scale > 1 (see Fig. A.1). All the quantities related to dust properties are computed using the optool¹ package (Dominik et al. 2021). Unbound grains sent on hyperbolic orbits because of radiation pressure are included in the simulation.

As will be explained in Section 2.2, the lifetime of the dust particles depends on the optical depth of the disk. This is simply parametrized by a maximum geometric optical depth τ_max which will be used to normalize optical depth profiles to estimate the collisional lifetime of the particles. As discussed in Zuckerman & Song (2012), the optical depth can be related to the fractional luminosity of the disk f_disk = L_disk/L_*, as τ_max ~ 2rf_disk/∆r, where L_disk is the luminosity of the disk, r is the radius of the disk and ∆r its radial width. The last parameters related to the properties of the birth ring are its inclination and position angle (i and φ, respectively). These two parameters relate only to the viewing geometry for computing the images (Section 2.3), without affecting the physics.

2.1.2 Gas parameters

The gaseous disk is parametrized by a normal distribution for the surface density distribution Σ_g(r), centered at a_g with a standard deviation σ_g, and normalized to a total gas mass M_gas. The composition of the gas is described with the mean molecular weight μ, and the gas temperature varies as T_g(r) = T_g,0(r/a_g)^−1/2.

2.2 Determining the dust density distribution

Since the gas exerts a force on the dust grains that depends on their velocities, we derive the dust density distribution in an iterative manner, as described below.

2.2.1 Inventory of forces

At t = 0, n_part particles are created and they immediately experience the following three forces: stellar gravity, radiation pressure, and gas drag. The first two are simply accounted for by reducing the stellar mass M_* by a factor (1 - β) when computing the acceleration (Burns et al. 1979). For the gas drag, we followed Jankovic et al. (2026) and assumed that the gas is in vertical hydrostatic equilibrium, with the vertical scale height computed as $$H_{\mathrm{g}}=c_{\mathrm{s}}(r)/{\mathrm{\Omega }}_{\mathrm{k}}=\sqrt{\frac{k_{\mathrm{B}}{T}_{\mathrm{g}}(r)r^3}{G{M}_{\star }\mu m_{\mathrm{H}}}},$$(3)

where c_s(r) is the local sound speed, Ω_k the Keplerian angular velocity in the midplane, k_B the Boltzmann constant, and m_H the mass of the hydrogen atom. The velocity of the gas in the vertical and radial directions are null, but because the gas feels its own pressure, its azimuthal velocity will depart from the local Keplerian velocity v_k and can be computed as $v_{\mathrm{g}}=v_{\mathrm{k}}\sqrt{1-\eta }$, where η is derived as in Jankovic et al. (2026, also capped at unity at large separations): $$\eta (r)=\frac{c_{\mathrm{s}}^2(r)r}{G{M}_{\star }}(\frac{7}{4}+\frac{r^2-r{a}_{\mathrm{g}}}{{\sigma }_{\mathrm{g}}^2}).$$(4)

At each distance r, there is a value of β(s) for which the effect of radiation pressure and gas drag cancel each other out (β(s) = η(r)), so that the particle no longer drifts in the radial direction. This divides the β(s) plane into two different domains of inward or outward drift. Figure A.2 shows four examples of the β(s) stability regions, determined by the condition η(r) = β(s), for different combinations of T_g,0 and μ. As discussed in Takeuchi & Artymowicz (2001, see also Krivov et al. 2009; Jankovic et al. 2026; Weber et al. 2026), in the absence of collisions, if a dust particle has a β value larger than the local value of η it will drift outwards, until it reaches a location in the disk where β(s) = η(r).

The thermal velocity of the gas is computed as in Takeuchi & Artymowicz (2001): $$v_{\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{m}}=\frac{4}{3}\sqrt{\frac{8k_{\mathrm{B}}{T}_{\mathrm{g},0}}{\pi \mu m_{\mathrm{H}}}{\left(\frac{r}{a_{\mathrm{g}}}\right)}^{-1/2}},$$(5)

The gas drag force can then be computed as: $$F_{\mathrm{g}\mathrm{a}\mathrm{s}}=-\pi {\rho }_{\mathrm{g}}{s}^2\sqrt{v_{\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{m}}^2+\mathrm{\Delta }{v}^2}\mathrm{\Delta }v,$$(6)

where ρ_g is the local gas density and ∆v is the relative velocity between the dust particle and the gas, the latter of which depends on η.

2.2.2 Integrating the equation of motion and assessing convergence

To follow the motion of the dust particles, we used a fourth order Runge-Kutta integrator and the acceleration accounts for all the aforementioned forces. The timestep δt is fixed to one hundredth of the orbital period at a distance of a_d. The final dust density distribution is estimated through an iterative process of several consecutive runs. The purpose of these runs is to make sure that the grain-to-grain collisions are accounted for properly. If one were to compute just one simple iteration, the end result would highly depend on its duration. We therefore followed the approach outlined in Thébault (2012), which is also similar to that described in Krivov et al. (2009).

A first iteration was computed (“run 0”) for an arbitrary duration (e.g., 500 000 years), which yielded a first estimate of the radial optical depth profile. For this first run, there are no collisions, and particles are not destroyed. It is parametrized by a total number of integrations, and at every timestep we increased the 1D vertical optical depth, τ₀(r), by the cross-section, πs², of the particles (meaning it is not computed only at the end of the simulation, but rather updated after each timestep). Once the run is finished, the τ₀(r) profile is normalized so that the maximum is equal to τ_max (an input parameter of the simulation). This is illustrated in Fig. 1, for which τ_max = 5 × 10⁻². We note that while the code is 3D, since the problem is axisymmetric we used azimuthally averaged quantities whenever possible (e.g., τ(r)) to speed up calculations.

The following iterations (n) will use the estimate of τ_n−1(r) of the previous iteration (n - 1) to determine whether a particle is destroyed or not. For each particle, at each timestep, we computed a probability that a particle is destroyed as $$f_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}}=\frac{\pi \delta \mathrm{t}{\tau }_{n-1}(r)}{t_{\mathrm{o}\mathrm{r}\mathrm{b}\mathrm{i}\mathrm{t}}},$$(7)

similarly to Olofsson et al. (2022). If a particle is destroyed, it is removed from the simulation and a collision does not produce new particles. The optical depth profile τ_n(r) is then updated, at each timestep, with the cross-sections of all remaining particles, and this is repeated until the simulation finishes. A run with collisions will be considered as finished when 99∙99% of the initial n_part have been destroyed. The inclusion of unbound particles requires an additional condition for particles to be removed from the simulation. These grains are set on hyperbolic orbits and will never come back in the birth ring. Their chances of being destroyed by collisions are therefore very small. We set a maximum distance (100 times the reference radius a_d) beyond which we considered the particles to be removed. Overall, the total duration of the simulation for the n^th iteration will depend on the collisional lifetime of the particles, since we needed to ensure the vast majority of particles are destroyed one way or another. It may be shorter or longer than in previous iterations, since τ(r) may differ, especially for the first few iterations. Figure 1 shows the evolution of the optical depth profiles for consecutive runs 0 to 5. It is apparent that the simulation quickly converges to a stable optical depth profile, where runs 3, 4, and 5 are almost perfectly identical to each other. For comparison, a slope of r^−1.5 is also displayed, which is the expected behavior for the optical depth beyond the birth ring (Strubbe & Chiang 2006; Thébault & Wu 2008). At large stellocentric distances, beyond ~2 000 au, the profile flattens towards r⁻¹ because of the increasing contribution of unbound grains (see Strubbe & Chiang 2006; Thébault et al. 2023).

Fig. 1 Geometric optical depth as a function of the stellocentric distance. A reference slope of r−1.5 is shown with a dotted line. The solid lines of various colors and thicknesses show the evolution of the optical depth for successive iterations (see Section 2.2.2 for details). The location of the birth ring is marked by the gray-shaded area (ad ± σd).

2.3 Synthetic images

During the last iteration (run 5), the code will also compute a set of images, in scattered light (total intensity or linear polarimetry) and thermal emission. The images are parametrized by a number of pixels and a pixel size. At each timestep, for each particle, we first applied two rotations to the (x, y, z) coordinates: the first one to account for the inclination i of the disk, and the second for its position angle φ. If the considered particle falls in the field of view of the image, the corresponding pixel is updated as described in the following.

2.3.1 Scattered light images

For each particle, we first computed the scattering angle (θ) as the arc cosine of the dot product between the line of sight and the position vector of the particle. Its contribution is then: $$F_{\mathrm{s}\mathrm{c}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{d}}=\frac{F_{\star }}{4\pi r^2}\pi s^2{S}_{11}(s,\theta )Q_{\mathrm{s}\mathrm{c}\mathrm{a}}(s)$$(8)

where S₁₁ is the element of the Müller matrix for the total intensity phase function, computed with optool (for polarized intensity this would be replaced by S₁₂), r is the distance to the star, and F_* the stellar flux.

2.3.2 Thermal emission

One key aspect when computing images in thermal emission is to estimate the temperature T_dust of the particles, which depends on the absorption coefficients (hence s and wavelength λ) and the distance r to the star. Assuming that the dust grains are at equilibrium temperature, we can equate how much light they absorb and emit, to obtain the following relationship: $$r(s,T_{\mathrm{d}\mathrm{u}\mathrm{s}\mathrm{t}})=\frac{R_{\star }}{2}\sqrt{\frac{{\int }_{\lambda }{F}_{\star }(\lambda )Q_{\mathrm{a}\mathrm{b}\mathrm{s}}(\lambda ,s)\mathrm{d}\lambda }{{\int }_{\lambda }\pi B(\lambda ,T_{\mathrm{d}\mathrm{u}\mathrm{s}\mathrm{t}})Q_{\mathrm{a}\mathrm{b}\mathrm{s}}(\lambda ,s)\mathrm{d}\lambda }},$$(9)

where B is the Planck function at the temperature T_dust and R_* the stellar radius. When initializing the simulation, the code samples an array of temperatures between 2 and 3000 K and computes the corresponding radius for all grain sizes (illustrated in Figure A.3). During the simulation, when computing the thermal emission image, for a particle of a given size and given distance, we inverted the above equation to find the temperature of the dust grain. If the particle falls in the field of view of the image, the contribution to the flux at a wavelength λ then is: $$F_{\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{l}}=4\pi s^2{Q}_{\mathrm{a}\mathrm{b}\mathrm{s}}(s,\lambda )\pi B(\lambda ,T_{\mathrm{d}\mathrm{u}\mathrm{s}\mathrm{t}}).$$(10)

2.3.3 Additional diagnostics

To facilitate the comparison between scattered light and thermal emission images, we can also save intermediate products, since we have the cross-section and distance for every particle in the simulation. For instance, we saved surface brightness profiles as a function of the radius for both kinds of images. This was done by summing the contributions of Equation (8) or (10) in radial bins (with Δr = 0.8 au) and dividing by the surface area of the considered ring. This is done regardless of the azimuth angles (and, hence, not at a fixed scattering angle in scattered light); this approach is similar to de-projecting the image and azimuthally averaging the surface brightness. The contributions to the radial profiles can also be saved for different intervals of β so that we can identify which grain sizes contribute most to the final image.

2.3.4 Convolution, sensitivity, and observables

For the purpose of this study, neither the images nor the radial profiles were convolved by a point spread function. Since we are mostly interested in the position of the maximum surface brightness, we also opted not to include noise in the final products. While these two aspects could not be overlooked if we were attempting to fit models to the observations, this is not the main goal of this work. Convolution would add further complexity to the multidimensional parameter space; for instance, even convolving the scattered light image with a point spread function of fixed width (e.g., 40 mas) would have different impacts depending on the distance of the star (~10 pc for a disk like AU Mic or 117.9 pc for HD121617).

Regarding the observables, some of the ARKS papers show radial profiles in surface density (e.g., Milli et al. 2026) instead of surface brightness. In this work, we decided that the final products should be as close as possible to the observations and therefore work with surface brightness profiles.

2.4 Caveats and known limitations

To be able to explore the parameter space, and explain a possible origin of the radial offsets, we needed to compute a rather large number of models. To maintain a reasonably short computing time, assumptions and simplifications have to be made. We mention here the main known limitations of the approach.

2.4.1 Collisional model

The collisional lifetime of the dust particles is a simplified estimation. For instance, the production of dust grains is assumed to follow a power-law distribution, and all subsequent collisions are purely destructive (the particle will simply be removed from the simulation).

Furthermore, in Equation (7), the relative cross-sections of the particles do not intervene when computing the probability that a particle is destroyed. In reality, we know that, for a s^−3.5 size distribution, the global geometrical cross section decreases as s^−0.5. As a consequence the collisional lifetime should increase as √s, thus making larger particles survive longer than smaller ones; this means that our collisional lifetime is underestimated. Additionally, according to Equation (7), the relative velocities between two particles are not accounted for when checking whether or not a particle is destroyed. The distribution of relative velocities between dust grains might play a more crucial role in simulations with gas compared to simulations without gas. On the one hand, gas drag should circularize the orbit (decreasing eccentricity while at the same time increasing pericenter distance, see Olofsson et al. 2022). This circularization is expected to lower the relative velocity between two particles since the velocity vectors should become more aligned with each other as the eccentricities of the particles decrease. On the other hand, most collisions take place in the birth ring (where the optical depth is maximum), where the velocities of the particles are the highest, regardless of their size. Given the complexity of implementing these two effects (cross-sections and velocities) and the underlying assumptions that would go behind such a model, we opted for a simpler solution that still accounts for the dynamical collision lifetime of the particles.

2.4.2 Gaseous disk

The prescription for the gaseous disk is also a simplistic description. We assumed a normal distribution for the radial distribution of gas. Such a distribution is qualitatively consistent with the observed radial distributions of CO emission in the ARKS sample (Mac Manamon et al. 2026) in the sense that CO density likely decreases both radially inwards and radially outwards from the belt center. We also assumed that the gas temperature is the blackbody one and is vertically constant. Furthermore, the gaseous disk has no velocity in the vertical direction, it does not viscously spread for the full duration of the simulation, and does not contribute gravitationally, despite being assigned a nonzero total mass. Marino et al. (2022) show that vertical mixing caused by turbulent diffusion could affect the spatial distribution of the dust particles, but since we have very few constraints on the mixing strength we opted not to include such effects in the simulation. Dust radial diffusion due to gas turbulence is not accounted for in the simulations as its strength remains unconstrained.

3 A first model

To have “realistic” default parameters for a first simulation, we chose HD 121617, an A-type star with a relatively narrow belt, as a proxy for a typical system hosting a gas-bearing debris disk. Following Marino et al. (2026), the stellar parameters are set to M_* = 1.9 M_⊙, L_* = 14 L_⊙, at a distance of 117.89 pc. The inclination and position angle of the disk are i = 44.1°, φ = 58.7°. Since we are not aiming to reproduce the ALMA and SPHERE observations of this system, the remaining parameters are defined more loosely. The birth ring of planetesimals is centered at a_d = 75 au with a standard deviation of σ_d = 5 au. The opening angle of the disk is set to ψ = 0.04 and the parent bodies have a proper eccentricity drawn from a normal distribution with a standard deviation of 0.05 (as in Thébault et al. 2023 for instance). The peak optical depth is set to τ_max = 5 × 10⁻³, corresponding to f_disk = 5 × 10⁻⁴ (assuming Δr = 3σ_d). The minimum and maximum grain sizes are 0.1 μm² and 5 mm, with a slope in q = −3.5, and β_scale = 1. The optical properties are computed with the distribution of hollow spheres model (DHS, Min et al. 2005) assuming a maximum filling factor (f_max) of 0.8 and the DIANA composition (Woitke et al. 2016 and references therein) for which ρ = 2.08 g cm⁻³. For these parameters, the size for which β = 0.5 is ~20 μm, meaning that a significant fraction of the size distribution (0.1 < s < 20 μm) consists of unbound grains. For the gaseous disk, the total mass is M_gas = 10⁻² M_⊕, with a distribution centered at a_g = a_d = 75 au, a standard deviation σ_g = 10 au, a mean molecular weight μ = 28 (implying a secondary origin for the gas, where the disk is CO rich and hydrogen poor), and a reference temperature of T_g,0 = 40 K (Brennan et al. 2026).

Figure 2 shows the collisional lifetime of particles as a function of the β value (only bound grains are considered in this map). When computing the last run of a simulation (run 5 from Fig. 1), at each timestep, if a particle with a given β value is destroyed, we incremented the corresponding (t,β) cell by unity. The color-coding in Fig. 2 therefore shows the number of collisions at a given (t,β), renormalized to the maximum for each β value. In this first simulation, most of the particles with β ≤ 0.1 will only survive for about 10 000 years. As mentioned in Section 2.4.1, the probability of having a destructive collision does not depend on the cross-section of the dust particles (Eq. (7)). The combined effect of the radiation pressure and gas drag is the only factor that could lead to a size-dependent collisional lifetime. Since particles with low β values are less sensitive to these forces, this yields a constant collisional lifetime. But, as β increases, so does the lifetime of the particles, since they spend a significant amount of time outside the birth ring, where densities are lower.

For this first model, we computed two images, one in scattered light, total intensity, at a wavelength of 1.63 μm and a second one in thermal emission at 880 μm³. As a benchmark, we also computed another simulation where the gas mass was set to 0. Figure 3 shows the results for the set of parameters described above, for the simulation with gas.

We can see that the radial profiles for the simulation with gas peak at the following distances: 78.5 au for SPHERE, and 75.5 au for ALMA (3.0 au offset). The results from the gas-free simulation show that the ALMA profile remains largely unaffected by gas drag, for this set of parameters. The SPHERE profile of the simulation with gas peaks ~1.7% further out relative to the simulation without gas (77.2 au), indicating that gas-dust interactions are shifting the distribution of small grains outward.

Figure 4 shows the cumulative contributions to the surface brightness as a function of the β value so that we can investigate which grain sizes are contributing to these radial profiles. This is a cumulative plot in the sense that the line β < 0.15 also includes the contribution of grains with β < 0.1 and so on. The lower panel is for the ALMA image and shows that the majority of the flux comes from grains with β < 0.1 (i.e., s > 100 μm) with a marginal contribution from grains smaller than that. In scattered light, the top panel of Fig. 4, the situation is quite different as all intervals contribute to the image. As β increases, the peak of the profile shifts to larger distances due to the combined effect of radiation pressure and gas drag. Interestingly, even though they are short-lived, unbound grains do contribute a non-negligible fraction to the image (they make up for the difference between the lines marked “All ß” and “ß < 0.5”), in line with the findings of Thebault & Kral (2019).

Fig. 2 Collisional lifetime as a function of the particles’ β value. The color-coding indicates the number of particles destroyed in a given cell, renormalized for each column.

Fig. 3 Results for the fiducial model, with a gas mass of 10−2 M⊕. Left : synthetic scattered light image in total intensity at 1.63 μm. Middle : thermal emission image at 880 μm. Right : normalized surface brightness profiles. The profiles for the simulation with (without) gas are shown with solid (dashed) lines. The ALMA profiles are shown in black and gray; the other two profiles (purple and orange) are for the SPHERE profiles.

4 Parameter space exploration

Even though our approach has relatively few free parameters, computing one simulation can be lengthy (from a few hours to a couple of days). We therefore could not compute a fully populated grid of models and could only explore combinations of a few parameters. We here followed an empirical approach to identify which ones drive the extent of the radial offset, guided by which ones can affect either the dynamics of the grains or the integrated optical properties. Table 1 summarizes the different families of models we computed, and the parameters that are varied with respect to the fiducial model are highlighted in bold font.

4.1 Gas mass and dust optical depth

The first parameters that we explored to vary the radial offset between the SPHERE and ALMA images were the gas mass (M_gas) and peak optical depth (τ_max) of the dust, since they will affect either the radial drift efficiency or the collisional timescale. Starting from the first model described in Section 3, we changed the gas mass M_gas to be between 10⁻³ and 1 M_⊕ (7 values). We also varied τ_max between 5 × 10⁻⁴ and 5 × 10⁻² (5 values), which resulted in 35 models (the other parameters remained the same as before; see family F1 in Table 1). For this set of parameters, the fractional luminosities of the disks would range between 5 × 10⁻⁵ and 5 × 10⁻³.

The peak positions and the offsets for these 35 models are displayed in Figure 5. In each panel, the five curves are for the different values of τ_max. For the low gas mass models, the peak positions of the SPHERE and ALMA profiles lie close to the reference radius a_d = 75 au. For both wavelengths, as the gas mass increases the peaks shift to larger stellocentric distances. We note that for the ALMA profiles (third panel) there is some dispersion on the peak positions (≲0.5 au). This is because the ALMA surface brightness profiles mostly trace the largest grains (Fig. 4), which are less numerous in our simulations given the steep size distribution, leading to larger numerical Poisson noise for the ALMA profiles compared to the SPHERE profiles.

As shown in the bottom panel of Figure 5, the offset is always positive (r_SPHERE ≥ r_ALMA). It first increases with the gas mass but tends to decrease as M_gas becomes larger, as even the large grains start to drift outward. For the higher gas masses, the offset increases significantly, reaching values close to 30 au. This happens when a “secondary ring” begins to appear in scattered light. Such secondary rings are known to arise from gas-dust interactions (Takeuchi & Artymowicz 2001; Jankovic et al. 2026). In brief, as discussed in Jankovic et al. (2026), the drift velocity of the dust particles, under the effects of stellar gravity, radiation pressure, gas drag, and Poynting-Robertson (PR) drag, strongly depends on their sizes. The location of the secondary ring is mostly determined by where the drift velocity of the particles falls to zero, which depends on gas pressure profile and density, which in turn depends on gas mass. In this work, we also accounted for collisions, and the cross-sections of the dust grains drive the shape of the optical depth radial profile, which in turn governs the collisional lifetime of the particles and, hence, the location of any secondary ring that may form, meaning that the location of the secondary ring also depends on τ_max. Since a simulation is considered as finished when the optical depth profile, τ(r), no longer changes between consecutive iterations, the secondary ring will not continue to evolve over time. Its brightness and detectability will, however, strongly depend on the crosssections of the grains, their scattering efficiencies (for scattered light images), and their absorption efficiencies and temperatures (for thermal emission images). Since this is investigated in depth in Jankovic et al. (2026), in the rest of this paper our main focus lies on disks with a single ring, but we still discuss our simulations that show secondary rings and the differences of our approach with the one of Jankovic et al. (2026) in Section 5.3.

Increasing gas mass and optical depth seems to have an opposite effect on the peak position of the radial profiles, which can be seen in the third panel of Figure 5. The curve for models with high optical depth (τ_max = 5 × 10⁻² in yellow) is below the one for models with low optical depth (τ_max = 5 × 10⁻⁴ in black). Regardless of the gas mass, the peak positions of the ALMA profiles peak farther out for low optical depth compared to larger optical depth; the dust particles have more time to drift before they are eventually destroyed by collisions. Furthermore, in the lower panel of Fig. 5, the curves for simulations with τ_max = 5 × 10⁻⁴ and 10⁻³ show a maximum offset for a gas mass of ~ 10⁻² M_⊕, and the offset decreases for higher gas masses (barring the last two points). As τ_max increases, the maximum offset is reached for higher gas masses, closer to 0.1 M_⊕, the offset then decreasing for higher gas masses. This correlation between the two parameters is further discussed in Section 5.1.

Fig. 4 Cumulative contributions to the surface brightness radial profiles as a function of β for SPHERE (top) and ALMA (bottom), for the fiducial model, with Mgas = 10−2 M⊕. The hatched area corresponds to ad ± σd.

Fig. 5 Peak positions for the radial profiles for SPHERE (first and second from the top) and ALMA (third from the top) and the difference between the two (bottom, in log-scale), as a function of the total gas mass and for different peak optical depths (legend in the middle panel), for family F1. The second panel from the top shows a zoomed-in version of the topmost panel. For the top three panels, the horizontal dashed lines show ad = 75 au.

4.2 The size distribution

The surface brightness in the SPHERE images strongly depends on the detailed dust properties, especially for the small end of the size distribution (Fig. 4). Since this would increase the parameter space significantly, we opted not to change the composition or porosity of the dust particles. To test if we could obtain larger offsets between SPHERE and ALMA, we varied parameters that are related to the particle size distribution, which will affect the Q_sca, S₁₁, and the cross-sections s². The first two parameters that we considered here are the slope of the size distribution q and the strength of the radiation pressure via β_scale.

First, changing the slope of the size distribution q (to values smaller than −3.5) will increase the number of small particles in the simulation. Such steeper values are indeed expected when taking into account the size-dependence of the crushing energy (Gáspár et al. 2012). Since the small particles are the ones most strongly affected by radiation pressure and gas drag, increasing their numbers relative to larger grains should affect the peak position of the scattered light surface brightness profile. Second, as mentioned before, the motivation to introduce the β_scale parameter is that it allows us to more finely control the β(s) relationship. From Eq. (1), this means that β(s) is no longer self-consistent with the dust properties Q_abs or Q_sca (or the stellar parameters L_* or M_*). In principle, changing β(s) to have a different value for s_blow can be achieved by changing either the dust composition or the porosity (e.g., Arnold et al. 2019), but since the Q_sca and S₁₁ strongly depend on the composition, this would not allow us to isolate the impact that s_blow has on the final images. Instead, by setting β_scale to ¼, the blow-out size decreases from 20 μm to 5 μm (F2 and F4), while all the other optical properties remain the same. This break of consistency with respect to Eq. (1) was deemed acceptable to limit the number of free parameters.

There is a third parameter, related to the size distribution, that we can change: the minimum grain size s_min (usually s_min < s_blow and all the grains with sizes in between would be gravitationally unbound). Family F5 is based on the parameters of F4, only with s_min changed from 0.1 μm to 0.5 μm.

We considered five families of models (F1 to F5). For each family, we varied M_gas and τ_max, as described in Table 1, but with q = −3.5 or −3.75 (close to the slope reported in Gáspár et al. 2012), β_scale = 1 or ¼, and with s_min = 0.1 or 0.5 μm. For each model we computed the offset between the peak positions of the SPHERE and ALMA profile. Then, for each family of models (F1 to F5), we computed a kernel density estimation (KDE), summing normal distributions with a standard deviation of 1 au, centered at the value of the offset. This is done to avoid sensitivity to the bin size and provide a more continuous representation of the data. The results are shown in Figure 6 for all nine families of models.

For some of the families, there are offsets larger than 15 au. This happens when the secondary ring becomes brighter than the birth ring in scattered light, and is further discussed in Section 5.3. Barring these particular cases, we see that the maximum of the KDEs shifts towards larger offsets as we decrease q and β_scale. Decreasing s_blow seems to be slightly more efficient than decreasing q; the KDE for family F2 peaks at larger distance than the one for family F3. Changing both parameters seems to be cumulative and the KDE for family F4 peaks at ~6.5 au. Releasing a larger quantity of small dust particles (steeper size distribution or decreasing the blow-out size via β_scale) therefore helps increase the offset between the SPHERE and ALMA radial profiles. We remark that increasing the minimum grain size from 0.1 to 0.5 μm (F5 compared to F4) yields a similar distribution of offsets. The scattering efficiencies Q_sca for submicron particles (s < 0.5 μm) are likely too small to have a significant impact on the resulting surface brightness profiles (see e.g., Figure 4 of Olofsson et al. 2023).

4.3 Gas kinetic temperature and mean molecular weight

Since the local sound speed c_s intervenes when estimating the gas-drag force and depends on both the temperature profile and the mean molecular weight, we then varied these two parameters. We computed new models based on family F4 (Table 1), which shows the largest offsets and set T_g,0 = 300 K and μ = 28 (F6) or T_g,0 = 40 K with μ = 2 (F7).

Figure 6 (for families F4, F6, and F7) shows that changing either parameter has little impact on the radial offset between SPHERE and ALMA. As illustrated in Fig. A.2, for a given radius r, the stable region where η(r) = β(s) (Eq. (4)) is reached for larger β values when increasing T_g,0 or decreasing μ. This means that, assuming they survive long enough, the small particles will not drift as far away from the birth ring. For the larger grains, contributing to the ALMA profiles, the situation is less straightforward. In Figure A.2, the stability curve for family F4 crosses r = 75 au for β ~ 10⁻³, which is shifted to β ~ 10⁻² for models of family F6. All particles with 10⁻³ < β < 10⁻² should drift outward for models of family F4, while the same grains should remain close to the birth ring for families F6 and F7. While this should help in obtaining a larger offset for families F6 or F7, the assumption that the particles can survive long enough may not hold. Indeed, the birth ring is where the optical depth peaks, hence where the lifetime of any particle is shortest. Overall, this is a multidimensional problem, but Fig. 6 suggests that the kinetic temperature or mean molecular weight have little impact on the offset between SPHERE and ALMA observations, in the regime r_SPHERE - r_ALMA ≲ 10 au. However, when it comes to the formation of secondary rings, Fig. 6 shows that increasing T_g,0 or decreasing μ results in fewer systems with offsets larger than ~20 au (that is, secondary rings). This is in line with the results presented in Jankovic et al. (2026) who found that the distribution of small particles is smoother for lower μ values, while grains with sizes 10 < s < 30 μm pile-up outside the birth ring and create a detectable secondary ring for models with larger μ values (their Figure 7).

Fig. 6 KDE and histograms of the offsets for all families of models (F1 to F9). The vertical dashed line is centered at 0. The fact that some of the curves display negative offsets is due to the fixed kernel’s width of 1 au.

4.4 Birth ring and gaseous disk radius

In this subsection, we investigate how the scale of the offset varies when changing the reference radius of either the birth ring a_d or the gaseous disk a_g. Our reference family of models is F4 (Table 1), since this family shows the largest offsets (Fig. 6). For family F8, the gaseous disk peaks at a_g = 70 au instead of 75 au, with a_d = 75 au (we note that for this family of models, μ is set to 2, compared to 28 for F4, but we show that this only has a small effect, Section 4.3). Family F9 is similar to F4, except that a_d = a_g = 50 au. Figure 7 shows the KDEs for these three families of models (with standard deviation of 0.05 for the kernel), but in this specific case, the offset is divided by a_g so that we can quantify the relative effect that a_g has. Even though the extent of the parameter space exploration is limited to a handful of values, our results suggest that the relative offset does not depend significantly on a_g since all KDEs peak at the same location. The offset would be smaller for more compact disks, and larger for disks at a greater separation from the star.

Fig. 7 KDE of the offset divided by the reference radius of the gaseous disk ag for families F4, F8, and F9, using a kernel standard deviation of 0.05.

5 Discussion

5.1 The competition between gas mass and dust optical depth

Figure 8 shows the offset measured between the peak positions of the SPHERE and ALMA radial profiles as a function of the ratio M_gas/τ_max. Since τ_max, which governs the collisional timescale, is directly related to the total dust mass, the ratio M_gas/τ_max can be considered as a proxy for the gas-to-dust mass ratio. We used Equation (7) from Wyatt (2008) to convert the fractional luminosity f_disk to an estimated dust mass, as $$M_{\mathrm{d}\mathrm{u}\mathrm{s}\mathrm{t}}=12.6f_{\mathrm{d}\mathrm{i}\mathrm{s}\mathrm{k}}{a}_{\mathrm{d}}^2{\kappa }_{\nu }^{-1}{X}_{\lambda }^{-1},$$(11)

where κ_ν is the dust opacity (1.9 cm² g⁻¹ at 0.89 mm, Marino et al. 2026, converted to 50.7 au² M_⊕⁻¹) and X_λ is meant to account for the emission falloff at longer wavelengths, which we assumed to be 4 as in Wyatt (2008). Since the linear transformation from τ_max to M_dust depends on a_d, we did not include results from family F9, since a_d is 50 au compared to 75 au for the rest of the models. The two populations display a similar pattern as a function of M_gas/τ_max, but with a vertical offset; models with smaller β_scale show overall larger offsets compared to models with β_scale = 1, as discussed in Section 4.2. We also show a 2D KDE underneath the scatter plot, for all the models, to better show the density distribution of offsets.

Initially, as the ratio M_gas/τ_max increases, so does the size of the offset. Higher gas masses will increase the offset as the gas drag force is larger (Equation (6)). Gas drag will both dampen the particles’ eccentricities and increase their pericenter distances faster for higher gas masses (e.g., Fig. 5 in Olofsson et al. 2022). They should thus be less likely to return inside the birth ring where their chances of being destroyed are higher. On the other hand, increasing τ_max will decrease the offset, since it influences the collisional lifetimes of the particles (Eq. (7)). For larger τ_max values, particles will be destroyed more efficiently, before they might reach the location where β = η. Lowering the maximum optical depth or increasing the gas mass therefore helps shift the peaks of the surface brightness profiles to larger separations.

However, beyond M_gas/τ_max ~ 10 (M_gas/M_dust ~ 0.3) the distribution of offsets splits into two different branches. The upper branch is populated with models that show a secondary ring brighter than the birth ring in scattered light; models with large gas masses or low optical depths where the small dust particles can drift outwards before being destroyed. On the other hand, the lower branch is populated with models where the effect of gas drag results in a broadening of the SPHERE surface brightness profiles, or the creation of a secondary ring that does not become brighter than the birth ring. As an example, the surface brightness profiles for a model with a “shoulder” in scattered light are shown in Figure 9. The contribution of the small grains in the SPHERE profile is mostly confined to distances in the range 90-100 au. As a consequence, the surface brightness in the birth ring (the hatched region) arises from larger grains, which are less sensitive to gas drag and radiation pressure, resulting in a smaller offset between SPHERE and ALMA.

Fig. 8 Offset between SPHERE and ALMA as a function of the ratio between the gas mass Mgas and maximum optical depth τmax for all models from families F1 to F8 (i.e., all but F9). Models with different βscale are shown with different symbols and colors. The dashed blue box shows the locus of models whose secondary belt is brighter than the birth ring in the SPHERE radial profiles. The underlying contours show a 2D KDE for all the models. A secondary x-axis at the top shows the estimated gas to dust mass ratio.

Fig. 9 Same as Figure 4 for the model of family F1 with Mgas = 10−2 M⊕ and τmax = 5 × 10−4.

5.2 Mid-IR observations

Mid-IR observations are complementary to mm observations as the thermal emission no longer probes the Rayleigh-Jeans tail and the flux density depends much more on the stellocentric distance. For each model we therefore computed new thermal emission images at a wavelength of 18 μm (akin to JWST/MIRI observations and in the same wavelength regime as ELT/METIS observations). The peak positions of the SPHERE, ALMA, and the newly computed radial profiles (to which we refer to as “JWST”) are shown in Figure 10 arranged by their respective families. The three rightmost panels show the histograms of peak positions for each wavelength.

In most cases, the peak positions of the JWST profiles lie in between the ALMA and SPHERE ones, but there are some exceptions. Interestingly, we see that the JWST profiles peak inward of the ALMA profiles for family F1. This can be explained by computing the flux density at 18 and 880 μm, which is shown in Figure A.4, for a given grain size (100 μm). The flux density decreases much faster with the stellocentric distance in the mid-IR compared to submillimeter. For the remaining families of models, this effect makes a smaller difference. The different behavior of family F1 compared to all the other models is related to the abundance of small dust grains. Family F1 is the only one for which q = −3.5 and β_scale = 1, and it is the only family for which the JWST profiles mostly peak inward of the ALMA profiles.

We can also see in Figure 10 that secondary rings can become brighter than the birth ring at mid-IR wavelengths, for families F6 to F8. The common denominator between these families are the gas properties, either the temperature (T₀) or the mean molecular weight (μ). The brightness ratio between any secondary ring and the main ring at mid-IR wavelength could therefore provide additional constraints on the origin of the gas, informing whether it is of primordial or secondary origin.

To better understand why the mid-IR peak positions seem to depend on the size distribution (e.g., families F1 vs. F4), Figure 11 shows the cumulative surface brightness profiles at 18 μm. Similarly to Fig. 4, the contributions of the different intervals of β are shown. For this example, most of the surface brightness can be explained by the contributions of either large grains (β < 0.05) or much smaller, unbound grains. As we varied either β_scale or q (or both; family F4) to have a larger number of small dust particles, the contribution of unbound grains increased, shifting the peak to larger separations. This is in line with the results presented in Thebault & Kral (2019), showing that the unbound grains can significantly contribute to the mid-IR flux. This underlines the importance of obtaining high angular resolution observations at these wavelengths, as they probe a population of dust particles that are otherwise poorly constrained.

Any further analysis of the mid-IR profiles would require the implementation of PR drag in the code. Recent JWST observations of debris disks such as Vega (Su et al. 2024) and Fomalhaut (Sommer et al. 2025) suggest that PR drag plays a significant role in shaping the morphology of the disks as seen at these wavelengths (even though these two particular systems have smaller fractional luminosities than the models discussed in this study).

5.3 Secondary rings

Milli et al. (2026) compared the radial profiles derived from SPHERE and ALMA observations for the ARKS sample and showed that the disk around HD 131835 is a rather unique system. Among the six disks that display a significant offset between scattered light and submillimeter observations, HD 131835 has the largest one: the ALMA profile peaks at 65 au (from the frank modeling results) to be compared to 112 au in scattered light for the brightest ring, leading to a relative offset of ~70%. Jankovic et al. (2026) investigated different scenarios to explain the origin of this large offset. One of them is the effect of gas drag, on which we based our calculations of η(r). However, there are some differences between the two approaches: Jankovic et al. (2026) did not directly account for collisions nor for the presence of unbound grains, while we did not include the effect of PR drag in our simulations.

Some of our models show first order similarities with the observations of HD 131835. Figure 12 shows one of the models from family F4, with a gas mass of 5 × 10⁻² M_⊕ and a maximum optical depth of τ_max = 5 × 10⁻⁴. This combination of parameters ensures that the grains can survive for a long time while experiencing non-negligible gas drag. The SPHERE, JWST, and ALMA images are displayed, from left to right, while the rightmost panel shows their corresponding surface brightness profiles. At mm wavelengths, only the main ring is visible, even though from the surface brightness there might be a hint of a secondary bump close to 100-110 au. At mid-IR wavelengths, the main ring appears very similar to the one seen with ALMA, but the secondary ring is almost as bright (~70%) in thermal emission and is clearly visible in the image. In scattered light, the secondary ring is brighter than the main ring, echoing what is observed for HD 131835 (Feldt et al. 2017; Milli et al. 2026; Jankovic et al. 2026). Interestingly, we can see in the image that the azimuthal brightness due to the scattering phase function is different in the inner and outer rings. This is because the inner ring mostly consists of large grains (β < 0.05) while the scattered light from the outer ring arises from particles with β > 0.05.

We note that HD 131835 is not the only debris disk for which a secondary ring has been reported from near-IR scattered light observations. Olofsson et al. (2023) reported the detection of a faint outer ring in the disk around HD 129590. This ring is only detected in total intensity, not in polarized light, and the disk is known to harbor cold CO gas (Kral et al. 2020). Similar results were found for HD 120326, with the exception that the disk seems to be devoid of cold CO gas (Bonnefoy et al. 2017; Desgrange et al. 2025). While it is often considered as a “hybrid” disk, the disk around HD 141569 is known to be gas-rich (Flaherty et al. 2016) while showing multiple rings in scattered light (Singh et al. 2021). However, none of these systems currently have high angular resolution observations at millimeter wavelengths that can be compared with the scattered light observations.

Interestingly, when the gas mass is large and the optical depth small, we find that some models also exhibit a secondary ring in the ALMA surface brightness profiles (though they never become brighter than the birth ring). Regardless of the wavelength at which a secondary ring is observed, such bimodal radial profiles could be misinterpreted as the presence of a gap in the distribution of planetesimals, mimicking the presence of massive perturbers in the disk. That being said, in a second-generation gas scenario, it is unclear how we could have, at the same time, a low optical depth and a high gas mass since the gas should be released from collisions of icy bodies; this would also produce significant amounts of dust (e.g., Marino et al. 2020).

These results should, however, be taken with caution since our simulations do not include the effect of PR drag. As mentioned in Jankovic et al. (2026), PR drag starts becoming important in regions where the gas density is low (but has little impact in regions of higher density). This is best illustrated in their Figure 7, showing grains can stop drifting outwards before they reach the regions where β(r) = η(r). In our simulations drift stops either at β(r) = η(r) or where the gas density and drift velocity drop to very low values. In practice, this suggests that the secondary ring could be closer to the birth ring compared to our results. Future works focusing on the properties of these secondary rings should include PR drag (as in Jankovic et al. 2026) and collisions (as in this study).

Fig. 10 Left : peak positions of the surface brightness for SPHERE, JWST, and ALMA for the different families of models listed in Table 1. The horizontal dashed line marks ad (75 au for families F1 to F8 and 50 au for family F9). Rightmost panels : histograms of the peak positions for each wavelength regime.

Fig. 11 Same as Figure 4, for the thermal emission surface brightness profile at 18 μm, for a model of family F1 (Mgas = 0.5 M⊕, τmax = 10−3).

6 Summary

With new high-angular-resolution observations we can probe the spatial distribution of dust particles at unprecedented details, opening new ways to study debris disks by following the dynamics of the dust grains. Motivated by the comparison between new ARKS submillimeter and archival SPHERE near-IR observations of gas-bearing debris disks (Milli et al. 2026), we investigated how the distribution of dust particles is affected by gas drag. The metric of interest in this study is the radial offset between peak surface brightness in images at different wavelengths. We have presented a new parameter space study, building upon past works (Takeuchi & Artymowicz 2001; Krivov et al. 2009; Olofsson et al. 2022; Jankovic et al. 2026), and we summarize our main results here:

• Gas drag can affect the dynamics of dust particles, and the final dust density distribution strongly depends on both the gas mass and the disk’s optical depth. Increasing the gas mass increases the rate at which particles drift outwards, while a greater optical depth implies shorter collisional lifetimes, limiting how far particles can drift.

• Even though we investigated only a small part of the parameter space, we find that variations in the gas mean molecular weight and kinetic temperature do not appear to have a significant impact on the extent of the offset;

• We find that either having a steeper size distribution (q < -3.5) or decreasing the blow-out size can increase the radial offset measured from multiwavelength observations, making it compatible with the offsets reported in Milli et al. (2026);

• We show that mid-IR wavelengths are complementary to near-IR and submillimeter observations and that measuring radial offsets between the three different combinations can provide additional constraints on the gaseous disk’s properties and the grain size distribution;

• While the contribution of unbound grains is negligible at submillimeter wavelengths, they do contribute in near-IR scattered light and even more so at mid-IR wavelengths, in thermal emission. Their contribution should not be ignored in this wavelength regime;

• When no secondary ring is observed, the extent of the offset scales with the reference radius of the disk: smaller disks display smaller offsets;

• Similarly to previous work we find that a secondary ring can appear outside the planetesimal belt. This only occurs under the right conditions (low optical depth and high gas mass). The ring appears in both scattered light and mid-IR observations (and is brighter than the birth ring), while it remains barely detectable at millimeter wavelengths. Models with a low mean molecular weight or a high kinetic temperature seem to be less prone to developing a secondary ring.

The approach described in this work is currently not well suited to directly fitting observations, due to the relatively high computational time required, nor regions of the parameter space that we have not explored (e.g., dust composition). That being said, the observations obtained with the ARKS Large Program, combined with near-IR scattered light observations, have opened a new avenue to studying the dynamics of dust particles, and this new diagnostic can help us better understand how dust particles evolve after being produced from collisions of planetesimals in debris disks.

Fig. 12 Left to right : SPHERE, JWST, and ALMA images, as well as surface brightness profiles for all three wavelengths, for a model from family F4, with a gas mass of 5 × 10−2 M⊕ and a maximum optical depth of 5 × 10−4. The pixel size is 12.26 milliarcseconds for all three images (akin to SPHERE/IRDIS), and no convolution by a point spread function has been applied.

Acknowledgements

We thank the referee for their thorough report that helped improving the presentation of our results, especially for the comparison between the different families of models. This research made use of Astropy (http://www.astropy.org) a community-developed core Python package for Astronomy (Astropy Collaboration 2013, 2018), Numpy (Harris et al. 2020), Mat-plotlib (Hunter 2007), and Numba (Lam et al. 2015). This paper makes use of the following ALMA data: ADS/JAO.ALMA# 2022.1.0033 8.L, 2012.1.00142.S, 2012.1.00198.S, 2015.1.01260.S, 2016.1.00104.S, 2016.1.00195.S, 2016.1.00907.S, 2017.1.00167.S, 2017.1.00825.S, 2018.1.01222.S and 2019.1.00189.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. The project leading to this publication has received support from ORP, that is funded by the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101004719 [ORP]. We are grateful for the help of the UK node of the European ARC in answering our questions and producing calibrated measurement sets. This research used the Canadian Advanced Network For Astronomy Research (CANFAR) operated in partnership by the Canadian Astronomy Data Centre and The Digital Research Alliance of Canada with support from the National Research Council of Canada the Canadian Space Agency, CANARIE and the Canadian Foundation for Innovation. MRJ acknowledges funding provided by the Institute of Physics Belgrade, through the grant by the Ministry of Science, Technological Development, and Innovations of the Republic of Serbia. SM acknowledges funding by the Royal Society through a Royal Society University Research Fellowship (URF-R1-221669) and the European Union through the FEED ERC project (grant number 101162711). MB acknowledges funding from the Agence Nationale de la Recherche through the DDISK project (grant No. ANR-21-CE31-0015). AMH acknowledges support from the National Science Foundation under Grant No. AST-2307920. SMM acknowledges funding by the European Union through the E-BEANS ERC project (grant number 100117693), and by the Irish research Council (IRC) under grant number IRCLA-2022-3788. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. EM acknowledges support from the NASA CT Space Grant. LM acknowledges funding by the European Union through the E-BEANS ERC project (grant number 100117693), and by the Irish research Council (IRC) under grant number IRCLA-2022-3788. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. JM acknowledges funding from the Agence Nationale de la Recherche through the DDISK project (grant No. ANR-21-CE31-0015) and from the PNP (French National Planetology Program) through the EPOPEE project. A.A.S. is supported by the Heising-Simons Foundation through a51 Pegasi b Fellowship. Support for BZ was provided by The Brinson Foundation. CdB acknowledges support from the Spanish Ministerio de Ciencia, Innovación y Universidades (MICIU) and the European Regional Development Fund (ERDF) under reference PID2023-153342NB-I00/10.13039/501100011033, from the Beatriz Galindo Senior Fellowship BG22/00166 funded by the MICIU, and the support from the Universidad de La Laguna (ULL) and the Consejería de Economía, Conocimiento y Empleo of the Gobierno de Canarias. JBL acknowledges the Smithsonian Institute for funding via a Submillimeter Array (SMA) Fellowship, and the North American ALMA Science Center (NAASC) for funding via an ALMA Ambassadorship. TDP is supported by a UKRI Stephen Hawking Fellowship and a Warwick Prize Fellowship, the latter made possible by a generous philanthropic donation. SP acknowledges support from FONDE-CYT Regular 1231663 and ANID - Millennium Science Initiative Program -Center Code NCN2024_001.

References

Appendix A Dust temperature, β, and η

We show here additional diagnostic plots providing context for the modeling approach. The β(s) and η(r) curves are shown in Figs. A.1 and A.2, respectively. Figure A.3 shows how the temperature of the dust particles varies as a function of the stel-locentric distance (r) and grain size (s) and Figure A.4 shows how the thermal emission flux density varies as a function of r for a specific grain size (s = 100 μm) at two different wavelengths (18 and 880 μm, respectively).

Fig. A.1 Radiation pressure strength β(s) for the default composition and stellar parameters (solid black line). Scaling down the original curve by βscale effectively changes the blow-out size sblow (red curve). The horizontal dashed line shows β = 0.5, while the vertical dashed lines show the blow-out size for both β(s).

Fig. A.2 Stability curves in the β(r) plane, for different configurations of the reference gas temperature and mean molecular weight. A particle released at a given (r,β) above (below) one of the curves will drift inwards (outwards), until it reaches the curve. The horizontal dashed line represents ad, the location of the birth ring.

Fig. A.3 Example map of Tdust as a function of the grain size and stel-locentric distance, with iso-temperature contours.

Fig. A.4 Flux density as a function of stellocentric distance, for a particle of size s = 100 μm (β ~ 0.1) for two different wavelengths, 18 and 880 μm.

All Tables

All Figures

	Fig. 1 Geometric optical depth as a function of the stellocentric distance. A reference slope of r−1.5 is shown with a dotted line. The solid lines of various colors and thicknesses show the evolution of the optical depth for successive iterations (see Section 2.2.2 for details). The location of the birth ring is marked by the gray-shaded area (ad ± σd).
In the text

	Fig. 2 Collisional lifetime as a function of the particles’ β value. The color-coding indicates the number of particles destroyed in a given cell, renormalized for each column.
In the text

Fig. 3 Results for the fiducial model, with a gas mass of 10−2 M⊕. Left : synthetic scattered light image in total intensity at 1.63 μm. Middle : thermal emission image at 880 μm. Right : normalized surface brightness profiles. The profiles for the simulation with (without) gas are shown with solid (dashed) lines. The ALMA profiles are shown in black and gray; the other two profiles (purple and orange) are for the SPHERE profiles.

In the text

	Fig. 4 Cumulative contributions to the surface brightness radial profiles as a function of β for SPHERE (top) and ALMA (bottom), for the fiducial model, with Mgas = 10−2 M⊕. The hatched area corresponds to ad ± σd.
In the text

Fig. 5 Peak positions for the radial profiles for SPHERE (first and second from the top) and ALMA (third from the top) and the difference between the two (bottom, in log-scale), as a function of the total gas mass and for different peak optical depths (legend in the middle panel), for family F1. The second panel from the top shows a zoomed-in version of the topmost panel. For the top three panels, the horizontal dashed lines show ad = 75 au.

In the text

	Fig. 6 KDE and histograms of the offsets for all families of models (F1 to F9). The vertical dashed line is centered at 0. The fact that some of the curves display negative offsets is due to the fixed kernel’s width of 1 au.
In the text

	Fig. 7 KDE of the offset divided by the reference radius of the gaseous disk ag for families F4, F8, and F9, using a kernel standard deviation of 0.05.
In the text

Fig. 8 Offset between SPHERE and ALMA as a function of the ratio between the gas mass Mgas and maximum optical depth τmax for all models from families F1 to F8 (i.e., all but F9). Models with different βscale are shown with different symbols and colors. The dashed blue box shows the locus of models whose secondary belt is brighter than the birth ring in the SPHERE radial profiles. The underlying contours show a 2D KDE for all the models. A secondary x-axis at the top shows the estimated gas to dust mass ratio.

In the text

	Fig. 9 Same as Figure 4 for the model of family F1 with Mgas = 10−2 M⊕ and τmax = 5 × 10−4.
In the text

	Fig. 10 Left : peak positions of the surface brightness for SPHERE, JWST, and ALMA for the different families of models listed in Table 1. The horizontal dashed line marks ad (75 au for families F1 to F8 and 50 au for family F9). Rightmost panels : histograms of the peak positions for each wavelength regime.
In the text

	Fig. 11 Same as Figure 4, for the thermal emission surface brightness profile at 18 μm, for a model of family F1 (Mgas = 0.5 M⊕, τmax = 10−3).
In the text

	Fig. 12 Left to right : SPHERE, JWST, and ALMA images, as well as surface brightness profiles for all three wavelengths, for a model from family F4, with a gas mass of 5 × 10−2 M⊕ and a maximum optical depth of 5 × 10−4. The pixel size is 12.26 milliarcseconds for all three images (akin to SPHERE/IRDIS), and no convolution by a point spread function has been applied.
In the text

	Fig. A.1 Radiation pressure strength β(s) for the default composition and stellar parameters (solid black line). Scaling down the original curve by βscale effectively changes the blow-out size sblow (red curve). The horizontal dashed line shows β = 0.5, while the vertical dashed lines show the blow-out size for both β(s).
In the text

	Fig. A.2 Stability curves in the β(r) plane, for different configurations of the reference gas temperature and mean molecular weight. A particle released at a given (r,β) above (below) one of the curves will drift inwards (outwards), until it reaches the curve. The horizontal dashed line represents ad, the location of the birth ring.
In the text

	Fig. A.3 Example map of Tdust as a function of the grain size and stel-locentric distance, with iso-temperature contours.
In the text

	Fig. A.4 Flux density as a function of stellocentric distance, for a particle of size s = 100 μm (β ~ 0.1) for two different wavelengths, 18 and 880 μm.
In the text