© The Authors 2025

1 Introduction

It has been well established that protoplanetary disks are composed of gas and a small fraction of dust, with a mass of about 1% of the total mass of the gas (e.g., Armitage 2024). However, a clear picture of the dynamics of these constituents has not yet emerged. For instance, there is still no general consensus on the possible mechanism that gives rise to the redistribution of angular momentum in protoplanetary disks. Over the last century, numerous efforts have been made to solve this problem. Currently, two mechanisms are considered as possible candidates to explain the transport of angular momentum in protoplanetary disks (see Turner et al. 2014; Lesur et al. 2023, for a review): (i) turbulence parameterized by a dimensionless α-parameter (Shakura & Sunyaev 1973), which can be driven by different sources (e.g., magnetorotational instability (MRI; Balbus & Hawley 1991), hydrodynamic instabilities such as the Rossby wave instability (RWI; Lovelace et al. 1999; Li et al. 2000, 2001), or the vertical shear instability (VSI; Arlt & Urpin 2004)); or (ii) magnetized disk winds (MDW; Blandford & Payne 1982).

In the frame of angular momentum transport driven by turbulence, the dynamics of the gas and dust in a protoplanetary disk are implicitly linked by the value of the α-parameter, since a non-zero value of α generates a certain level of turbulence in the gas which in turn induces turbulent diffusion (i.e., stirring) of the dust through the action of drag forces (Youdin & Lithwick 2007). Nevertheless, it is important to mention that dust coupling to gas turbulence can lead to diffusion (Chen & Lin 2020) and anti-diffusion (i.e., dust clumping, see e.g., Johansen et al. 2007). Both processes can affect the formation of planets (e.g., Venturini et al. 2020). Indeed, turbulent dust diffusion can prevent dust accumulation in regions of a protoplanetary disk where density or pressure gradients are at work.

With all the above in mind, we are interested in studying the effects of turbulent dust diffusion on planetary migration. In recent 2D numerical studies, where dust diffusion was not taken into account, it has been proposed that nonsymmetric dusty structures may arise in the vicinity of a planet (Benítez-Llambay & Pessah 2018; Regály 2020; Guilera et al. 2023; Chrenko et al. 2024). The so-called dust-void and dust filament structures (regions of low and high dust density, respectively) can produce a total torque on the planet that in many cases (depending on the Stokes number, St, of the dust and on the dust-to-gas mass ratio, ϵ), is positive and therefore leads to an outward migration for low-mass planets (0.3 ≤ M_p ≤ 10 M_⊕, with M_⊕ the Earth mass).

In this work, we examine the effect of turbulent dust diffusion on the density distribution of dust around a low-mass planet through multifluid simulations. In particular, we wish to analyse the formation of the dust-void and filament under the action of the turbulent diffusion of dust. We began with 2D models, similar to those considered by Benítez-Llambay & Pessah (2018), to quantify the change in the torque exerted on the planet when dust can diffuse. Subsequently, we investigate the possible implications of 3D dust distributions on planetary migration. We did so by performing 3D multifluid simulations, parameterizing the dust sedimentation through the δ-diffusion coefficient (Morfill & Voelk 1984; Weber et al. 2019). Overall, this study is aimed at determining whether the formation of a 3D dust-void and a filament structure are possible and whether gravitational torques are affected.

The paper is laid out as follows. In Sect. 2, we describe the physical model of a dust-gas disk. The code and numerical setups used in our 2D and 3D multifluid simulations are provided in Sect. 3. In Sect. 4, we show the results of our 2D numerical models on the dust-void destruction due to turbulent dust diffusion. The dust-void morphology obtained from 3D simulations and the effects on the migration of a low-mass planet are presented in Sects. 5 and 6. Finally, a brief discussion and our conclusions are given in Sect. 7.

2 Physical model

In this section, we describe different components of our physical model, including turbulent dust diffusion and the back-reaction of dust on gas employed in this study.

2.1 Governing equations

We considered a 3D, non-self-gravitating gas-dust disk whose evolution is governed by the following equations: $${\mathrm{\partial }}_t{\rho }_{\mathrm{g}}+\mathrm{\nabla }\cdot \left({\rho }_{\mathrm{g}}\mathbf{v}\right)=0,$$(1) $${\mathrm{\partial }}_t\left({\rho }_{\mathrm{g}}\mathbf{v}\right)+\mathrm{\nabla }\cdot ({\rho }_{\mathrm{g}}\mathbf{v}\otimes \mathbf{v}+p\mathbf{l})=-{\rho }_{\mathrm{g}}\mathrm{\nabla }\mathrm{\Phi }+\mathrm{\nabla }\cdot \tau -{\mathbf{f}}_{\mathrm{d}},$$(2) $${\mathrm{\partial }}_t{\rho }_{\mathrm{d}}+\mathrm{\nabla }\cdot ({\rho }_d\mathbf{u}+{\mathbf{j}}_{\mathrm{d}})=0,$$(3) $${\mathrm{\partial }}_t\left({\rho }_{\mathrm{d}}\mathbf{u}\right)+\mathrm{\nabla }\cdot ({\rho }_{\mathrm{d}}\mathbf{u}\otimes \mathbf{u})=-{\rho }_{\mathrm{d}}\mathrm{\nabla }\mathrm{\Phi }+{\mathbf{f}}_{\mathrm{d}},$$(4)

where ρ_g, ρ_d, v, and u denote the gas density, dust density, gas velocity, and the dust velocity, respectively. Furthermore, Φ denotes the gravitational potential, | is the unit tensor, and p is the gas pressure. For the latter, we consider the locally isothermal equation of state: $$p=c_s^2{\rho }_{\mathrm{g}},$$(5)

where c_s is the isothermal sound speed. Note that the sound speed is kept constant, namely, we did not include an energy equation. Consequently, the aspect ratio of the gas disk, h_g, remains constant. The aspect ratio of the dust, h_d, is initially set as a constant (Eq. (22)), at t = 0.

To mimic the turbulence in the gaseous component, we included a viscous stress tensor in the right-hand side (RHS) of the Eq. (2), whose 3D form can be expressed as $$\tau ={\rho }_{\mathrm{g}}\nu [\mathrm{\nabla }\mathbf{v}+(\mathrm{\nabla }\mathbf{v}{)}^{\mathrm{T}}-\frac{2}{3}(\mathrm{\nabla }\cdot \mathbf{v})\mathrm{l}],$$(6)

where v is the gas kinematic viscosity, parameterized through the usual Shakura & Sunyaev (1973) description $$\nu =\alpha c_s{H}_{\mathrm{g}}.$$(7)

In the equation above, H_g is the pressure scale height of the gas in the disk, defined as H_g = c_s/Ω_Kep where Ω_Kep is the Keplerian angular frequency. It is important to note that the α-parameter, apart from playing a role in the dynamics of the gas, also affects the dynamics of the dust because it generates a constraint on dust diffusion, as we describe below.

2.2 Dust diffusion flux

To take into account the possible effect of gas turbulence on the dust dynamics, in the dust mass conservation equation (see Eq. (3)), we include the gradient in the dust concentration through the vector j_d, given as (for details, see Morfill & Voelk 1984): $${\mathbf{j}}_{\mathrm{d}}=-D_{\mathrm{d}}({\rho }_{\mathrm{g}}+{\rho }_{\mathrm{d}})\mathrm{\nabla }\left(\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}+{\rho }_{\mathrm{d}}}\right).$$(8)

The quantity D_d is the dust diffusion coefficient, which is related to gas turbulence through the Schmidt number $$\mathrm{S}\mathrm{c}=\frac{\nu }{D_{\mathrm{d}}},$$(9)

which quantifies the relative effectiveness of the gas angular momentum transport and the dust-mixing processes. In this study, we consider Sc = 1, since we are analyzing a dust fluid with a Stokes number below unity (see Youdin & Lithwick 2007). Note that although both diffusive coefficients (v and D_d) have a common origin in the turbulence in the disk (and in this particular case, they are even equal), they describe different physical processes. Therefore, in order to differentiate these two quantities, we chose to use the dimensionless diffusivity parameter $$\delta =\frac{D_{\mathrm{d}}}{c_s{H}_{\mathrm{g}}},$$(10)

which is analogous to the dimensionless α-parameter.

2.3 Dust back-reaction force

In Eqs. (2) and (4), the last term on the RHS is the drag force term given as $${\mathbf{f}}_{\mathrm{d}}=\frac{{\mathrm{\Omega }}_{\mathrm{K}\mathrm{e}\mathrm{p}}}{\mathrm{S}\mathrm{t}}{\rho }_{\mathrm{d}}(\mathbf{v}-\mathbf{u}),$$(11)

where St is a dimensionless measure of the dust-gas coupling that relates to the dust grain properties via $$\mathrm{S}\mathrm{t}=\sqrt{\frac{\pi }{8}}\frac{r_{\mathrm{m}}{\rho }_{\mathrm{m}}\mathrm{\Omega }}{{\rho }_{\mathrm{g}}{c}_s},$$(12)

where r_m is the grain size and ρ_m is the material density of the dust. Note that the effect of the feedback force term given in Eq. (11) becomes more relevant in 3D disk models, where dust settling leads to large concentrations of dust at the midplane of the disk.

3 Multifluid hydrodynamical simulations

To numerically solve Eqs. (1)–(4), we used the publicly available hydrodynamic code FARGO3D¹ (Benítez-Llambay & Masset 2016; Benítez-Llambay et al. 2019) with an implementation of the fast orbital advection algorithm (FARGO; Masset 2000) used in our 2D models, along with the new implementation of the rapid advection algorithm on arbitrary meshes (RAM; Benítez-Llambay et al. 2023) employed in our 3D simulations.

In both 2D and 3D multifluid simulations, the dust density profile is initialized from the gas density using the constant dust-to-gas mass ratio, ϵ, and composed of separate fluids with Stokes numbers of St = 0.04, 0.26 and St = 0.55 at r₀ = 1. These Stokes numbers fall into the gas-dominated, transitional, and gravity-dominated regimes, respectively (see cases a), b) and c) analyzed in Benítez-Llambay & Pessah 2018). Although for small values of the Stokes number and gas viscosity (hence, for small values of δ, such as 10⁻⁴ used here), the feedback may be ineffective, when St approaches unity the feedback can produce significant changes in the radial velocity of the gas (Kanagawa et al. 2017). Therefore, to get an idea of the effect of turbulent dust diffusion on the dust structures around of the planet, in case c) we did not include the feedback.

3.1 Setup of 2D simulations

Our 2D numerical simulations were carried out following the framework of Benítez-Llambay & Pessah (2018) (see their Sect. 2), which, in turn, is based on the gas-dust disk model of Weber et al. (2018); only with the two major differences in the disk model mentioned above: we included the effect of the back-reaction of the dust on the gas and the diffusion of the dust. The effects of these two processes, neglected in the study of Benítez-Llambay & Pessah (2018), are presented in Sect. 4. Below, we briefly recall some aspects of the initial conditions.

The 2D multifluid equations (i.e., 2D version of the set of Eqs. (1)–(4)) are solved on a polar grid centred on to the star and in a frame co-rotating with the planet of orbital radius, r_p. The surface density of the gas disk is chosen to be a power law of r: $${\mathrm{\Sigma }}_g(r)={\mathrm{\Sigma }}_0{\left(\frac{r}{r_p}\right)}^{-\sigma },$$(13)

with Σ₀ = 2 × 10⁻³/π in code units (which corresponds to ≈ 200gcm⁻² in cgs units at r_p = 5.2 au) and a surface density slope of σ = 1/2. The temperature decays as a power law of the radius, as well, with a slope of β = 1. We use h_g = 0.05 and consider α ∈ [0, 3 × 10⁻³]. The dust density profile is initialized from Σ_g, using the constant dust-to-gas mass ratio ϵ = Σ_d/Σ_g = 0.01.

The gravitational potential Φ is composed of the gravitational potential of the star, the planet, and the indirect term (arising from non-inertial forces, i.e., the second term on the RHS in Eq. (19)). Additionally, in the gravitational potential of the planet, we applied a softening length, r_sm, equal to the planet Hill radius, r_H = r_p(M_p/3M_⋆)^1/3 (here, M_p and M_⋆ are the planet and star masses, respectively) to avoid numerical issues. Note that the value of r_sm, for the low-mass planets considered in this work, is significantly smaller than the value required in a gaseous disk to produce an approximate agreement between 2D and 3D torques (e.g., Tanaka & Okada 2024); whereas it is comparable to or larger than the thickness of the dusty disk. This entails that we overestimate the gaseous torque and underestimate the dusty torque, so that our statements about the importance of the dusty torque are conservative. We applied the same cut-off of the torque inside the Hill sphere and we also considered the same smoothing length as in Benítez-Llambay & Pessah (2018) for comparison purposes, since our objective is to describe the effects of diffusion and feedback on dust distribution. However, we are aware that this dust distribution is sensitive to the choice of the smoothing length in 2D (Regály 2020; Guilera et al. 2023), whose value is not well-determined (for a detailed discussion see Chrenko et al. 2024; Hou & Yu 2025).

Unless otherwise specified, the planet has a mass of M_p = 1.5 M_⊕ and is held on a fixed circular orbit around a star of a mass of M_⋆ = 1 M_⊙. We modeled the entire ring of a disk (2π in azimuth) with a radial extent between r = 0.48 and r = 2.08, considering N_r = 3130 and N_ϕ = 12 288 radial and azimuthal uniform grid zones, respectively.

3.2 3D disk simulations

The setup of our 3D model is based on the gaseous disk of Masset & Benítez-Llambay (2016, as described in their Appendix A). We used spherical coordinates (r, θ, ϕ), where r is the radial distance from the star, θ is the polar angle ($\theta =\frac{\pi }{2}$ at the midplane of the disk), and ϕ is the azimuthal angle. The aspect ratio of the gas disk is h_g = H_g/r. The gas disk density is given by $${\rho }_{\mathrm{g}}={\rho }_{\mathrm{g}}^{\mathrm{e}\mathrm{q}}(\sin \theta {)}^{-\beta -\xi +h_{\mathrm{g}}^{-2}},$$(14)

with $${\rho }_{\mathrm{g}}^{\mathrm{e}\mathrm{q}}=\frac{{\mathrm{\Sigma }}_0}{\sqrt{2\pi }{h}_{\mathrm{g}}{r}_p}{\left(\frac{r}{r_p}\right)}^{-\xi },$$(15)

where Σ₀ is the surface density at r = r_p. We set β = 1, ξ = 1.5 and h_g = 0.05 in all models considered here. We initialize the gas velocity components as follows: v_r = v_θ = 0 and $$v_{\varphi }=\sqrt{\frac{G{M}_{\star }}{r\sin \theta }-\xi c_s^2}.$$(16)

The gravitational potential Φ is given by $$\mathrm{\Phi }={\mathrm{\Phi }}_S+{\mathrm{\Phi }}_p,$$(17)

where $${\mathrm{\Phi }}_S=-\frac{G{M}_{\star }}{r},$$(18)

and $${\mathrm{\Phi }}_p=-\frac{G{M}_p}{\sqrt{r^{\mathrm{\prime }2}+r_{\mathrm{s}\mathrm{m}}^2}}+\frac{G{M}_{p}r\cos (\varphi -{\varphi }_p)\sin \theta }{r_p^2}$$(19)

represent the stellar and planetary potentials, respectively. In Eq. (19), r′ = |r − r_p| is the cell-planet distance, with the planet located at $(r_p,{\theta }_p,{\varphi }_p)=(1,\frac{\pi }{2},0)$ and r_sm is softening length of the planet’s potential. The second term on the right-hand side of Eq. (19) is the indirect term arising from the reflex motion of the star. Our simulations were performed with r_sm = 0.05H_g, which is comparable to twice the cell size of our numerical grid (we carried out other experiments with a larger r_sm and found similar results). In our 3D models, no cut-off of the torque is applied inside the Hill sphere.

Since dust grains above the midplane are not pressure-supported, the vertical gravitational force drives their settling on a timescale that is much shorter than that of viscous gas evolution. Consequently, the scale heights of dust and gas can be different (h_d ≠ h_g). Therefore, we considered an initial dust density profile that depends on the dust pressure scale h_d, such that $${\rho }_{\mathrm{d}}={\rho }_{\mathrm{d}}^{\mathrm{e}\mathrm{q}}(\sin \theta {)}^{-\beta -\xi +h_{\mathrm{d}}^{-2}}$$(20)

and $${\rho }_{\mathrm{d}}^{\mathrm{e}\mathrm{q}}=ϵ\frac{{\mathrm{\Sigma }}_0}{\sqrt{2\pi }{h}_{\mathrm{d}}{r}_0}{\left(\frac{r}{r_p}\right)}^{-\xi }.$$(21)

In Eq. (21), the initial dust-to-gas mass ratio ϵ is such that Σ_d/Σ_g = 0.01. The dust pressure scale is defined through the dimensionless diffusivity parameter δ as $$h_{\mathrm{d}}=h_{\mathrm{g}}\sqrt{\frac{\delta }{\delta +\mathrm{S}\mathrm{t}}}.$$(22)

The radial and vertical components of the dust velocity are both zero initially, whereas the azimuthal component is set to start from a Keplerian rotation profile. Contrary to the 2D models, the planet is allowed to migrate in all of our 3D simulations. To include the aerodynamic back-reaction (feedback) from the dust on the gas (cases with St = 0.04 and St = 0.26), we performed simulations using two fluids at a time (gas plus one dust species).

3.2.1 Computational grid

Due to the computational cost of 3D multifluid simulations, to achieve a numerical resolution similar to that used in the 2D models (see Sect. 3.1), we used the Fargo3D code with the RAM module. This module allows us to use a nonuniform mesh, with smaller cells near the planet (Benítez-Llambay et al. 2023). This approach helps to increase the resolution around the dust-depleted region at a more affordable computational cost. The numerical grid covers the disk region [r_min, r_max] = [0.48r_p, 2.08r_p] in the radial directions, $[{\theta }_{min},{\theta }_{max}]=[\frac{\pi }{2}-h_p,\frac{\pi }{2}]$, in co-latitude (we only simulated one hemisphere of the disk and assumed symmetry in the other). The azimuthal extent is [ϕ_min, ϕ_max] = [−π, π], covering the whole extent of the disk.

Under this numerical scheme, the bumps in the radial and azimuthal mesh functions (see Appendix B in Benítez-Llambay et al. 2023) are parameterized by considering the constants a = 0.2618, b = 0.3141 to model the resolution transitions, and with c = 1.5 and c = 15 for the radial and azimuthal resolution contrast between the two levels, respectively. The vertical resolution and spacing were chosen on the basis of the dust aspect ratio (see Eq. (22)). When h_d ≥ 0.1h_g, we used a uniform grid with constant resolution in the polar direction. Setting the number of grid cells as (N_r, N_θ, N_ϕ) = (1864, 100, 1864) leads to cube-sized cells Δr = r_pΔθ = r_pΔϕ = 5 × 10⁻⁴ at the planet’s position. We also note that the grid attains a resolution of H/100. Otherwise, if h_d < 0.1h_g, we used a mesh density function in the polar direction, whose resolution ratio is specified in Appendix A.

3.2.2 Boundary conditions

We used boundary conditions for each dust species similar to those implemented for the density and velocity components of the gas (see Chametla & Masset 2021, for details). We used damping boundary zones, as in de Val-Borro et al. (2006), with the width of the inner and outer damping rings being 0.15r_p and 0.5r_p, respectively. The damping timescale at the edge of each ring is equal to 0.3 of the local orbital period. We emphasize that the radial extension of the domain is large enough to avoid any disturbances from the disk boundaries affecting the dust-void region. Any perturbations found in this region are therefore purely a result of the dynamic interaction between the gas, dust, and the planet. Additionally, we used damping zones only in the radial direction because their inclusion in the vertical direction would drastically modify dust settling and might lead to an overly large dust-to-gas mass ratio in the midplane. Lastly, since we consider only one hemisphere of the disk here, we used reflecting boundary conditions at the midplane.

4 Destruction of the 2D dust-void in a gaseous turbulent disk

In this section, we describe how the 2D dust density structures near a M_p = 1.5 M_⊕ planet change when the effects of turbulent diffusion and dust feedback are included. For the purpose of maintaining a clearer notation throughout this study, we omit stating the value of α when δ is specified because these parameters are identical.

4.1 Effect of turbulent dust diffusion

To gain insight into the effect of turbulent dust diffusion, we first present the results of the models with St = 0.55, where the dust feedback on the gas is neglected and the influence of the turbulent flow diffusion can be isolated more clearly.

In Fig. 1, we show the relative perturbation of dust density and the dust streamlines around the planet in the case of a Stokes number of St = 0.55 and six different values of δ, calculated at t = 50 orbits. When δ = 0, a dust-void structure appears behind the planet, this feature is encircled by a high-density ring-shaped region, while in front of the planet a filamentary structure is formed. Looking at the other panels, we find that the dust-void structure and the dust filament undergo minor modifications for values of δ < 3 × 10⁻⁴.

As the viscosity and turbulent dust diffusion increase, the distribution of dust around the planet is altered. When δ ≥ 3 × 10⁻⁴, the size of the dust-void is reduced and the filament in front of the planet starts to vanish. Given a diffusivity parameter value of δ = 3 × 10⁻³, we find that neither the dust-void structure nor the filament form. This last case deserves particular attention since Benítez-Llambay & Pessah (2018) reported the formation of the dust-void and filament for the gas viscosity α = 3 × 10⁻³ (and δ = 0; see their Figure 1).

This result is noteworthy because the value of α used by Benítez-Llambay & Pessah (2018) is actually comparable with recent estimates of turbulence levels inferred from observed disks (Flaherty et al. 2020; Jiang et al. 2024). Therefore, turbulence dust diffusion can be critically important in these situations. It should be emphasized that we recovered the same structure as case c) reported in Benítez-Llambay & Pessah (2018) when we considered α = δ = 0, that is, in a dusty gaseous disk without turbulence (see top-left panel of Fig. 1).

The trend appearing in Fig. 1 can be understood as a competition between two antagonic processes. The distribution of the dust results from its advection (an inward drift and a horseshoe U-turn) and diffusion due to its interaction with the turbulent gas. The dust settles into a pattern whose characteristics depend on the timescale for diffusion (which tends to fill the dust-void) and that of advection (which tends to result in a dust-void). When the timescale for the former is short enough, namely, when δ is sufficiently large, the dust-void is suppressed.

Fig. 1 Relative perturbation of dust density for a Stokes number St = 0.55 (representative of the gravity-dominated regime), generated by an Earth-mass planet (Mp = 1.5 M⊕). The solid gray lines show the dust streamlines. All images are shown for a time of t = 50 orbits. It should be noted that when the dust diffusion is not included (upper left panel), the formation of the dust-void behind the planet and the dust filament in front of the planet exhibit the same structure, as reported in Fig. 1 (case c), of Benítez-Llambay & Pessah (2018). However, when dust diffusion increases, the dust-void and filament eventually disappear (see bottom-right panel).

4.2 Lengthening and weakening of high-density dust-void barrier: Dust back-reaction force effect

Similarly, for a Stokes number St = 0.26 (transitional regime not shown here), we find again that the dust-void and the filament structure begin to disappear for δ ≥ 3 × 10⁻⁴ and are absent when δ ≃ 3 × 10⁻³. We emphasize that when the dust-void and filament structures are maintained, which happens for δ < 3 × 10⁻⁴, their shape and extension resemble those found in the gravity-dominated regime. In other words, the radius of the ring-shaped high-density structure that delimits the dust-void increases as the diffusivity parameter decreases. Dust from the surroundings cannot penetrate this region because of aerodynamic friction by the gas and scattering by the planet. Additionally, we find that the high-density filament in front of the planet becomes more prominent as the parameter δ is reduced.

In general, the explanation for the widening of the dust-void and the elongation of the filamentary structure lies in the fact that the feedback and the turbulent diffusion of the dust incidentally conspire to counteract the scattered dust flow generated by the planet. That is to say, turbulent diffusion produces a spreading of the ring-shaped structure, but due to the dispersion produced by the planet, the dust tends only to diffuse away from the planet. However, because the dust exerts a drag force on the gas, the pressure of the gas increases in the region where the dust has moved to, generating a new pressure maximum that prevents the dust located at the outer edge of the ring from escaping. This results in a ring-shaped structure with a larger radius. A similar effect occurs in the filament in front of the planet: while turbulent diffusion tries to spread the dust, the scattering generated by the planet forces the dust to move away from it, but the feedback counteracts it, leading to an elongation of the structure in question.

For the case when the dust is strongly coupled to the gas (case St = 0.04, i.e., gas-dominated regime), we find that the dust-void and filament structures are drastically modified compared to the case when turbulent dust diffusion and dust feedback are discarded (see first and second panels in Fig. 3). For the case of δ = 0, the dust barrier delimiting the dust-void is dragged toward the edges of the gas horseshoe region, where there may be vortices propagating because the gas disk is inviscid (since α = 0). A similar effect occurs in the dust filament formed in front of the planet, since some of the dust approaching the filament from upstream is carried toward the edge of the horseshoe region (see second panel in Fig. 3). Therefore, the dust-void is left unprotected, which leads to a shallower dust-void and the dust filament is weakened.

Lastly, if there is a non-zero turbulent diffusion of the dust, the blurring of the dust-void and the filament is stronger and they are completely erased when δ = 3 × 10⁻³ (see lower panel in Fig. 3).

Fig. 2 Temporal evolution of the dust (top row) and total (bottom row) torques on a planet of mass Mp = 1.5 M⊕ for representative cases of the gas-dominated (left), transitional (middle) and gravity-dominated (right) regimes.

4.3 Impact of the dust-void on the planet’s torque

An immediate consequence of the dust-void and filamentary structures not surviving turbulent dust diffusion is reflected in the torque that the dust exerts on the planet. In the case of a gravity-dominated regime (St = 0.55, without feedback), unlike Benítez-Llambay & Pessah (2018), we find that the dust torque on the planet is negative when α = 3 × 10⁻³. Note that the total torque is also negative (see the third column in Fig. 2). For a lower level of gas turbulence (which implies a lower value of δ), the dust torque becomes less negative and it is possible to obtain a null total torque when δ ≃ 3 × 10⁻⁴. Reducing the diffusivity further, we find that for δ < 3 × 10⁻⁴, the total torque on the planet is always positive. This latter result is relevant in regimes of low turbulence in the gas of protoplanetary disks (e.g., Pinte et al. 2016; Lesur et al. 2023; Jiang et al. 2024), as it entails that planetary migration can be directed outward.

In the case of the transitional regime (St = 0.26, with feedback included), due to the broadening and lengthening of the dust-void and filament, respectively, we find similar values of the dust torque as in the gravity-dominated regime. For instance, we refer to cases with δ = 10⁻⁴ and δ = 3 × 10⁻³ of the second and third columns in Fig. 2. However, it is worth mentioning that in the interval of 3 × 10⁻⁴ ≤ δ ≤ 10⁻³, the dust torque is negative for St = 0.26, which differs from the behavior of the torque in the gravity-dominated regime over the same range of δ values (where the dust torque is positive and its magnitude increases when δ → 0; see Fig. 2).

Finally, it is important to note that when the barrier of the dust-void is broken or destroyed, which occurs in the case of the gas-dominated regime (St = 0.04), the dust torque is positive only for very small values of the δ-parameter (~10⁻⁵). Therefore, even in the presence of fine dust, the total torque on the planet is negative.

5 3D dust-void morphology

Here, we focus on the results of the 3D multifluid simulations. In addition to the case of M_p = 1.5 M_⊕ considered in Sect. 4, we extended our analysis toward more massive planets. Additionally, we let the planets migrate instead of keeping them on fixed orbits.

5.1 Dust-void destruction

In Fig. 4, we show the relative perturbation of dust density around a planet of M_p = 1.5 M_⊕ at the midplane of the disk for the gas-dominated regime (St = 0.04). We have chosen three different values of the dust turbulent diffusion δ = 10⁻⁵, 10⁻⁴, and 3 × 10⁻³ as reference cases, allowing us to reveal the main differences in the dust density distribution, as compared to 2D models.

When δ = 10⁻⁵, the dust-void starts to show a shape resembling its 2D counterpart. However, the ring-shaped high-density structure at the edge of the dust-void is spread out, mainly due to the effect of turbulent dust diffusion. Note that even though the feedback effect is less efficient in this case because the dust is initially confined to a very thin layer (h_d ≪ h_g), there is a sufficient contribution from this to maintain dust accumulation, thereby leading to a process similar to that explained in the Sect. 4.2. Also, the ring-shaped high-density structure exhibits perturbations characteristic of an epicyclic motion when r > r_p (see Morbidelli & Nesvorny 2012). Additionally, we find that the filament in front of the planet is erased and a trench emerges in its place.

In the case of δ = 10⁻⁴, the dust-void forms again but, unlike the 2D dust-void (see Figs. 1 and 3), it exhibits low central density and is no longer delimited by a ring-shaped barrier of high density. Additionally, the trench and the epicyclic oscillations (in both structures) arise as for δ = 10⁻⁵. Note that the edges of the dust-void and trench are wider compared to those of the ring-shaped barrier and filament regions of the 2D models. This is likely a consequence of turbulent dust diffusion, which spreads out dust density gradients.

Lastly, when δ = 3 × 10⁻³, the size of the dust-void and the trench are considerably reduced, which means that both structures are filled by “fresh” dust. Seen from another perspective, dust scattering is only efficient very close to the planet, resulting in a negligible dust-void and shallower trench. The absence of the dust-void and filamentary structures for δ = 3 × 10⁻³ is similar in both 2D and 3D models (see bottom panels in Figs. 3 and 4).

Fig. 3 2D dust-void structure for three different values of δ in the gas-dominated regime (St = 0.04), for a planet with Mp = 1.5 M⊕. Top: case without dust turbulent diffusion in a viscous gas disk, which is similar to the case a) studied in Benítez-Llambay & Pessah (2018, their figure 1). Middle: case of an inviscid gaseous disk (α = 0) without dust turbulent diffusion (δ = 0). Bottom: dust-void and filament are destroyed for δ = 3 × 10−3. Snapshots at t = 50 orbits.

Fig. 4 3D dust-void structure for three different values of δ in the gas-dominated regime (St = 0.04), for a planet with Mp = 1.5M⊕. Top: snapshot (calculated at t = 70 orbits) in the case of δ = 10−5. A low-density trench forms in front of the planet instead of a filament. Middle: case with δ = 10−4. The dust-void structure resembles the shape found in 2D simulations but without the high-density ring-shaped barrier. Bottom: dust-void and filament disappearing for δ = 3 × 10−3 (snapshot calculated at t = 100 orbits).

Fig. 5 Relative perturbation of dust density at the midplane of the disk for α = δ = 10−4 (top) and α = δ = 3 × 10−3 (bottom) calculated at t = 100 orbits. Three different planetary masses in the transitional regime (St = 0.26) are considered (see the labels of individual panels). Note: the planet is located at the center of each plot.

5.2 Deformation and pollution of the dust-void for planets with a higher mass

Next, we analyzed the evolution of dust structures in 3D disks when the mass of the planet increases. For this proposal, we ran a set of simulations with M_p = 3, 6 and 12 M_⊕. We set the Stokes number to St = 0.26.

In Fig. 5, we show the results of these numerical experiments. The relative perturbation of dust density at the midplane for a planet with M_p = 3 M_⊕ exhibits a well-defined structure for the dust-void and the filament, for both values of the dust diffusion coefficients. However, a detailed inspection of the δ = 10⁻⁴ case shows that the filament becomes a trench that connects directly to the dust-void. Another interesting feature is the overdense stripes adjacent to the dust-void and the trench in front of it. This dust pattern may be a consequence of buoyancy resonances (Chametla et al., in prep.) since an isothermal dusty gaseous disk behaves like an adiabatic pure gaseous disk (Lin & Youdin 2017), in which the buoyancy resonances produce a similar stripe-like pattern in gas (see Ziampras et al. 2023, and references therein).

For a planet of M_p = 6 M_⊕, when δ = 10⁻⁴, the trench in front of the planet opens up enough to cover a considerable region compared to the dust-void, which in turn becomes polluted with material orbiting at r > r_p. Otherwise, the dust-void and the trench structures are well defined when δ = 3 × 10⁻³. That is to say, there is no material entering the dust-void region and the width of the trench does not increase as in the previous, lower mass case. Lastly, for the most massive planet considered in this study (i.e., M_p = 12 M_⊕) and a value of δ = 10⁻⁴, we find again that the trench in front of the planet widens and it also contains some filaments of dust within it. Interestingly, an overdense island forms within the dust-void, which survives until the end of the simulation. However, at a larger turbulent dust diffusion δ = 3 × 10⁻³, the dust-void remains clear and the trench is divided from the dust-void by a very thin filament that emerges from the planet position.

6 Runaway and oscillatory migration driven by the dust-void evolution

6.1 Runaway migration signature

The runaway migration of planets in a gas disk (without dust component) is a self-sustaining mechanism with short migration time-scales ≤O(10²) orbits (Masset & Papaloizou 2003) and it applies mainly to super-Earth, sub-Jovian, and Jupitertype planets (Masset & Papaloizou 2003; Pepliński et al. 2008a,b,c; Lin & Papaloizou 2010; Paardekooper 2014). The main engine of the gas-driven runaway migration is a non-linear corotation torque maintained by the change in the fluid elements’ orbital radius executing U-turns at the end of horseshoe orbits. This torque scales with the rate of migration of the planet and with the vorticity-weighted coorbital mass deficit, δm (Masset & Papaloizou 2003). In the high-mass regime, whose deficit essentially corresponds to a coorbital mass deficit, which is proportional to the difference between the surface density trapped in the coorbital region and that of the orbit-crossing flow, while in the low mass regime, it essentially amounts to a jump of vortensity between the horseshoe region and the surrounding disk (Paardekooper 2014; McNally et al. 2018).

In the first regime, two necessary conditions for runaway migration to occur in a gaseous disk are that the planet creates a partial gap and that the mass deficit δm ≈ M_p. However, note the masses of the planets (M_p ∈ [1.5, 12] M_⊕) considered in this study are not large enough to satisfy both conditions simultaneously. Therefore, the migration of these low-mass planets in a gas disk must fall within the type I migration regime (see Chametla et al. 2023, and references therein). As a proof of this last argument, in the upper panel of Fig. 6 we show the temporal evolution of the semimajor axis of a planet with M_p = 1.5 M_⊕ obtained from a 3D dust-free model. Clearly, the planet’s semimajor axis, a, does not change in a runaway fashion; rather, it varies at a constant rate over the first 150 orbital periods, as in the type I migration regime over small orbital displacements.

The situation can be much different when a planet migrates in a dusty gaseous disk. We first analyze the migration for the case of M_p = 1.5 M_⊕ and St = 0.55. In the bottom panel of Fig. 6, it can be seen that the semimajor axis of the planet changes considerably over a period of time that is less than a hundred orbital periods. The latter suggests that the planet is experiencing a runaway migration. To confirm this hypothesis, Fig. 7 shows the temporal evolution of the planet’s drift rate when the planet is embedded in a gas and gaseous-dusty disk with δ = 10⁻⁴ and δ = 3 × 10⁻³, respectively. Clearly, the drift rate of the planet migrating in a gaseous-dusty disk is higher and varies considerably over the first orbital periods, even if the semimajor axis itself does not vary significantly during that short time frame.

It should be noted that while the first points in Fig. 7 show a systematic increase in the migration velocity over time, characteristic of a runaway migration, it tends toward an asymptotic value over a time comparable to the libration time, τ_lib ≡ 4πa/(3/2)Ωx_s (with x_s half-width of horseshoe region); here, on the order of 60–70 orbits. The first phase corresponds to a runaway migration; whereas in the final stage, migration settles into a type III regime. Also note that the velocity at which migration saturates in type III is the relative velocity between the planet and the dust. Since dust (especially St = 0.55) migrates rapidly inward (u_r ~ −1 × 10⁻³), the planet’s inward migration must be faster for the feedback to be broken. This suggests that the drift rate is not a function of the semimajor axis only, but also of the drift rate itself, as occurs in runaway migration. Note that the drift rate in a gaseous disk is constant which is to be expected due to the behavior of the planet’s semimajor axis (see upper panel in Fig. 6).

6.2 Planetary drift rate as oscillatory-torque migration trace in dusty gaseous disks

We return to the analysis of the migration in the cases when St = 0.04 and St = 0.26, as shown in the second and third panels of Fig. 6 for a planet of M_p = 1.5 M_⊕. Except for the case of δ = 3 × 10⁻³ in the gas-dominated regime (St = 0.04), which exhibits a almost stagnant oscillating migration, we found that in these models the behavior of the semimajor axis is far from type I migration and in turn does not change as quickly while maintaining the same direction as in the case of runaway migration found in the gravity-dominated regime.

A detailed inspection of each of the semimajor axes shows a migration with several changes of direction in short time compared to the dynamical timescale, which resembles stochastic migration (Nelson & Papaloizou 2004; Nelson 2005; Baruteau & Lin 2010). However, it is important to emphasize that there is no global turbulence on the dust distribution. The changes in migration are mainly due to deflection and deformation of the dust-void and the trench, respectively, producing an oscillatory behavior in the torque felt by the planet.

Note that due to the time over which the 3D multifluid simulations are run, a convergence study of the running time-averaged torque is beyond the scope of this work. Therefore, we henceforth prefer to adopt the term “oscillatory-torque migration” to refer to such changes in planetary migration, rather than stochastic migration.

Remarkably, for planets with a larger mass (M_p ∈ {3, 6, 12} M_⊕) in the transitional-dominated regime (St = 0.26) also we find an oscillatory-torque migration. Figure 8 shows the temporal evolution of the semimajor axis, the migration rate, and the eccentricity for these models. In the first column of Fig. 8, it can be seen that the temporal evolution of the planet’s semimajor axis exhibits frequent oscillations within an orbital period for both values of dust turbulent diffusion. However, there appears to be a marked trend with respect to the mass of the planet, since for δ = 10⁻⁴, as the mass of the planet increases, migration goes from outward to inward migration. On the other hand, for δ = 3 × 10⁻³, even when the mass of the planet increases, the migration is directed outwards.

As a proof of oscillatory behavior in planet migration, in the second column of Fig. 8 we show the drift rate of these planets with a larger mass. In these plots we have included the sign of the drift rate, since if there are mixed positive and negative values over the dynamical time it means that the planet experiences episodes of oscillatory migration. In all cases, we find positive and negative values in the drift rate of the planet. However, particular attention must be paid to the case of M_p = 6 M_⊕ and δ = 3 × 10⁻³, we find that for t ≤ 70 orbital periods, the semimajor axis grows rapidly in a runaway fashion. We also observe that the sign of the drift rate does not change on this time interval. Therefore, this suggests that for this particular case, the planet first experiences a runaway migration and subsequently an oscillatory-torque migration.

Lastly, in the third column of Fig. 8, we show the temporal evolution of the eccentricity. It can be seen that although there is an oscillatory behavior, there is also a marked increasing tendency for the eccentricity in almost all models. Except in the case when M_p = 12 M_⊕ and δ = 10⁻⁴ where there is an inward oscillatory-torque migration scenario. Interestingly, these results suggest that an outward oscillatory-torque migration inevitably results in the development of planetary eccentricity.

Fig. 6 Temporal evolution of the semimajor axis of a planet with Mp = 1.5 M⊕ embedded in a purely gaseous 3D disk (top) and in 3D dusty disks representative of the gas-dominated, transitional, and gravity-dominated regimes ((second, third, and bottom panel, respectively).

Fig. 7 Temporal evolution of the migration rate $\dot{a}$ for a migrating planet of Mp = 1.5 M⊕ embedded in a gaseous disk and also in a dusty-gaseous disk with δ = 10−4 and δ = 3 × 10−3 for a gravity-dominated regime (St = 0.55) without any feedback of the dust onto the gas.

6.3 Triggering oscillatory-torque migration

As we describe above, a planet migrating in a disk of gas and dust can experience either a runaway migration or an oscillatorytorque migration. Both types of migration are governed by the level of turbulent diffusion and feedback of the dust which shapes the structure of the dust-void and the trench. For instance, in the gas-dominated regime (St = 0.04) when δ = 10⁻⁵, the dust-void is partially surrounded by broken ring-shaped high-density structure and the trench shows an elongated pattern (see Fig. 4). These features lead to an inward oscillatory-torque migration. On the other hand, for δ = 10⁻⁴, we find that the dust-void and the trench feature survive, driving an outward oscillatory-torque migration. Finally, for δ = 3 × 10⁻³, the dust-void and the trench are significantly weakened and the planet experiences a nearly stagnant oscillatory-torque migration.

Furthermore, in the case of St = 0.55 (gravity-dominated regime), where the back-reaction force of the dust on the gas is not included, we find that the dust-void and the trench maintain a well-defined shape, resulting in an inward runaway migration for both values of δ = 10⁻⁴ and δ = 3 × 10⁻³, respectively. Therefore, it appears from the models described above that the outcome of the runaway migration is mainly governed by the evolution of the dust-void. In fact, unless the dust-void is deformed or populated with dust islands, runaway migration occurs. Otherwise, any modification of the dust-void inevitably results in an oscillatory-torque migration.

7 Discussion and conclusions

In this work, we report the results of 2D and 3D multifluid hydrodynamics simulations, including dust turbulent diffusion, aerodynamic dust feedback, and planet migration (in the 3D case). Our aim was to determine how the inclusion of aforementioned effects changes the so-called dust-void and filament dust structures formed behind and in front of low-mass planets.

The importance of analyzing the possible effects of turbulent dust diffusion and feedback on these structures in the dust lies in the fact that Benítez-Llambay & Pessah (2018) recently showed that the asymmetric distribution of dust around low-mass planets embedded in 2D disks can produce a torque that dominates the other components of the gas torques (corotation and Lindblad torques). In a subsequent study, Guilera et al. (2023) quantified the impact due to the distribution of dust on the migration of the planets. However, in both studies, 2D laminar disks were considered, discarding the effect of turbulent diffusion of dust as well as the dust feedback. It is worth mentioning that in both studies the viscosity of the gas parameterized by the alpha coefficient was included, which was set to α = 3 × 10⁻³.

With this in mind, we started our study with 2D models for a low-mass planet (M_p = 1.5 M_⊕) on a fixed circular orbit embedded in a dusty-gaseous disk following the same framework as in Benítez-Llambay & Pessah (2018). Nevertheless, here we have taken into account the turbulent dust diffusion (parameterized by δ ∈ [0, 3 × 10⁻³]) and the back-reaction force from the dust on the gas. We find that the dust-void and the filament are formed when δ ≲ 10⁻⁴. On contrary, higher levels of turbulent dust diffusion prevent the formation of structures around the planet, reducing the ability of dust-driven torques to cause outward migration.

Of course, there could be other factors contributing to the modification of the 2D dust-void apart from the turbulent diffusion and feedback of dust; for instance, the smoothing radius of the planet’s potential (Regály 2020) or the absence of pebble accretion (Chrenko et al. 2024). The first possibility, beyond its role as a physical variable, arises as a limitation of 2D numerical models since it is chosen ad hoc. Therefore, in order to obtain a more appropriate physical description of the evolution of the dust structures, a 3D numerical treatment is necessary.

Consequently, we included in our study 3D simulations of migrating planets with M_p ∈ [1.5, 12] M_⊕ (which also includes the representative case analyzed in the 2D models). These simulations reveal that the morphology of the dust-void and the filament depends strongly on the turbulent dust diffusion as well as on the inclusion of feedback. The changes in the dust structures can be summarized as follows:

• ■

  In the case of smaller values of dust turbulent diffusion, δ < 10⁻⁴, the dust-void is drastically modified and the ring-shaped high-density structure surrounding it is weakened and scattered. Note that these structures still do not resemble the structures found at Benítez-Llambay & Pessah (2018), where neither turbulent dust diffusion nor feedback were included.

• ■

  For a δ value of the order of ~10⁻⁴, the dust-void survives and its fate will simply depend on how much dust falls into the dust-void.

• ■

  In the case of δ values larger than 10⁻⁴, the dust-void is almost filled for lowest mass planet and for larger masses, the dust-void decreases considerably.

• ■

  Instead, for all values of δ, the high-density filament in front of the planet is replaced by a trench.

Remarkably, we find runaway migration or an outward (inward) oscillatory-torque migration driven by the dust-void evolution due to the gas feedback and turbulent dust diffusion. It must be emphasized that the evolution of the dust-void and trench depends on the δ value, the dust feedback, and (slightly) on the mass of the planet. Lastly, the direction of the runaway migration (inward or outward) depends on the evolution of the trench, which is intrinsically linked to the turbulent diffusion and feedback of the dust.

Fig. 8 Temporal evolution of the orbital parameters for the planets of larger mass migrating in a 3D dusty disk including δ = 10−4 and δ = 3 × 10−3 for the transitional regime with St = 0.26. Left. Semimajor axis a. Center. Migration rate $\dot{a}$. Right. Planet’s eccentricity, e.

Acknowledgements

We thank the anonymous referee for his constructive and careful report, which significantly improved the quality of this manuscript. This work was supported by the Czech Science Foundation (grant 21-23067M). The work of O.C. was supported by the Charles University Research Centre program (No. UNCE/24/SCI/005). Computational resources were available thanks to the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID: 90254).

Appendix A Vertical dust height resolution treatment

From Eq. (22), it follows that the vertical height of the dust depends on the Stokes number and the coefficient of turbulent diffusion of the dust. Therefore, we have that for the pair of values of (St, h_d) = (0.04, 3 × 10⁻³) and (St, h_d) = (0.55, 10⁻⁴) studied in our dusty gaseous disks, the aspect ratio of the dust takes the respective maximum and minimum values $h_{\mathrm{d}}^{max}$ = 0.013 and $h_{\mathrm{d}}^{min}$ = 6.7 × 10⁻⁴. This last value of h_d is of the order of the resolution value with cube-sized cells Δr = r_pΔθ = r_pΔϕ = 5 × 10⁻⁴ at the planet position when a uniform mesh is used in the vertical component.

To avoid possible numerical effects on the dust dynamics around the planet, for the resolution of vertical component of the computational mesh we have adopted two prescriptions: 1) a uniform vertical domain with a constant resolution Δθ = 5 × 10⁻⁴ and, 2) A non-uniform vertical domain where the resolution increases considerably for values of h_d < 0.1h_g such that, it is satisfied that h_d/Δθ ≥ 10.68. In Fig. A.1 we show the resolution ratio at colatitude used in 2). Within the domain of the computational mesh between $\frac{\pi }{2}-0.01$ and $\frac{\pi }{2}$, the resolution is kept constant. To maintain the aforementioned constraint, the number of zones, N_θ, is modified according to the value of h_d in each model.

We note that due to the high computational cost of our 3D models when increasing the number of zones in the vertical component of the grid, we have fixed the h_d/Δθ ratio to a little more than ten grid cells for h_d. Nevertheless, we think that it is a good compromise to keep this minimum value for the numerical resolution of h_d. In Fig. A.2 we show the comparison of the temporal evolution of the semimajor axis of a planet with M_p = 3 M_⊕ for the pair of parameters (St, h_d) = (0.26, 10⁻⁴), applying both grid configurations. It is clear that the change is minimal between the two prescriptions of the vertical grid. Therefore, it is safe to adopt this vertical resolution ratio in description 2).

Fig. A.1 Colatitude resolution ratio for the models with a dust aspect ratio less than 0.1hg. Note: in the domain between $\theta =\frac{\pi }{2}-0.01$ and $\theta =\frac{\pi }{2}$, the resolution is constant and it satisfies the minimum requirement hd/Δθ = 10.86 in each model.

Fig. A.2 Comparison of the temporal evolution of the semimajor axis of a planet of mass Mp = 3.0M⊕ for α = δ = 10−4 and St = 0.26 using prescriptions 1) and 2) for the resolution of the grid in the vertical component (magenta and purple lines, respectively).

References

All Figures

Fig. 1 Relative perturbation of dust density for a Stokes number St = 0.55 (representative of the gravity-dominated regime), generated by an Earth-mass planet (Mp = 1.5 M⊕). The solid gray lines show the dust streamlines. All images are shown for a time of t = 50 orbits. It should be noted that when the dust diffusion is not included (upper left panel), the formation of the dust-void behind the planet and the dust filament in front of the planet exhibit the same structure, as reported in Fig. 1 (case c), of Benítez-Llambay & Pessah (2018). However, when dust diffusion increases, the dust-void and filament eventually disappear (see bottom-right panel).

In the text

	Fig. 2 Temporal evolution of the dust (top row) and total (bottom row) torques on a planet of mass Mp = 1.5 M⊕ for representative cases of the gas-dominated (left), transitional (middle) and gravity-dominated (right) regimes.
In the text

Fig. 3 2D dust-void structure for three different values of δ in the gas-dominated regime (St = 0.04), for a planet with Mp = 1.5 M⊕. Top: case without dust turbulent diffusion in a viscous gas disk, which is similar to the case a) studied in Benítez-Llambay & Pessah (2018, their figure 1). Middle: case of an inviscid gaseous disk (α = 0) without dust turbulent diffusion (δ = 0). Bottom: dust-void and filament are destroyed for δ = 3 × 10−3. Snapshots at t = 50 orbits.

In the text

Fig. 4 3D dust-void structure for three different values of δ in the gas-dominated regime (St = 0.04), for a planet with Mp = 1.5M⊕. Top: snapshot (calculated at t = 70 orbits) in the case of δ = 10−5. A low-density trench forms in front of the planet instead of a filament. Middle: case with δ = 10−4. The dust-void structure resembles the shape found in 2D simulations but without the high-density ring-shaped barrier. Bottom: dust-void and filament disappearing for δ = 3 × 10−3 (snapshot calculated at t = 100 orbits).

In the text

	Fig. 5 Relative perturbation of dust density at the midplane of the disk for α = δ = 10−4 (top) and α = δ = 3 × 10−3 (bottom) calculated at t = 100 orbits. Three different planetary masses in the transitional regime (St = 0.26) are considered (see the labels of individual panels). Note: the planet is located at the center of each plot.
In the text

	Fig. 6 Temporal evolution of the semimajor axis of a planet with Mp = 1.5 M⊕ embedded in a purely gaseous 3D disk (top) and in 3D dusty disks representative of the gas-dominated, transitional, and gravity-dominated regimes ((second, third, and bottom panel, respectively).
In the text

	Fig. 7 Temporal evolution of the migration rate $\dot{a}$ for a migrating planet of Mp = 1.5 M⊕ embedded in a gaseous disk and also in a dusty-gaseous disk with δ = 10−4 and δ = 3 × 10−3 for a gravity-dominated regime (St = 0.55) without any feedback of the dust onto the gas.
In the text

	Fig. 8 Temporal evolution of the orbital parameters for the planets of larger mass migrating in a 3D dusty disk including δ = 10−4 and δ = 3 × 10−3 for the transitional regime with St = 0.26. Left. Semimajor axis a. Center. Migration rate $\dot{a}$. Right. Planet’s eccentricity, e.
In the text

	Fig. A.1 Colatitude resolution ratio for the models with a dust aspect ratio less than 0.1hg. Note: in the domain between $\theta =\frac{\pi }{2}-0.01$ and $\theta =\frac{\pi }{2}$, the resolution is constant and it satisfies the minimum requirement hd/Δθ = 10.86 in each model.
In the text

	Fig. A.2 Comparison of the temporal evolution of the semimajor axis of a planet of mass Mp = 3.0M⊕ for α = δ = 10−4 and St = 0.26 using prescriptions 1) and 2) for the resolution of the grid in the vertical component (magenta and purple lines, respectively).
In the text