Quantifying the Impact of the Dust Torque on the Migration of Low-mass Planets

Disk solids are critical in many planet formation processes, however, their effect on planet migration remains largely unexplored. Here we assess for the first time this important issue by building on the systematic measurements of dust torques on an embedded planet by Benitez-Llambay&Pessah (2018). Adopting standard models for the gaseous disk and its solid content, we quantify the impact of the dust torque for a wide range of conditions describing the disk/planet system. We show that the total torque can be positive and revert inward planet migration for planetary cores with $M_{\rm p} \lesssim 10 M_\oplus$. We compute formation tracks for low-mass embryos for conditions usually invoked when modeling planet formation processes. Our most important conclusion is that dust torques can have a significant impact on the migration and formation history of planetary embryos. The most important implications of our findings are: $\it{i})$ For nominal dust-to-gas mass ratios $\epsilon \simeq 0.01$, low-mass planets migrate outwards beyond the water ice-line if most of the mass in solids is in particles with Stokes numbers St $\simeq 0.1$. $\it{ii})$. For $\epsilon \gtrsim 0.02-0.05$, solids with small Stokes numbers, St $\simeq 0.01$, can play a dominant role if most of the mass is in those particles. $\it{iii})$ Dust torques have the potential to enable low-mass planetary cores formed in the inner disk to migrate outwards and act as the seed for massive planets at distances of tens of au.


INTRODUCTION
As planets grow and evolve in the protoplanetary disks where they form, they modify the environment surrounding them, which in turn can alter the planet's evolution.In particular, a planet embedded in a disk leads to asymmetric structures with the ability to exert a net torque on the planet that can alter its dynamics.The ensuing planet migration is a key process in planet formation (see Venturini et al. 2020c, and references therein).While the presence of disk solids is widely recognized as critical in many processes involved in planet formation and evolution, its role has not been assessed in planet migration.Building on the first systematic measurements of dust torques on an embedded planet reported by Benítez-Llambay & Pessah (2018, hereafter BLP18), this paper is a first attempt to address this important issue.
In sufficiently viscous disks, the gas torque exerted onto the planet has two important contributions: the differential Lindblad torque and the corotation torque.The former is negative for typical protoplanetary disk models while the latter depends sensitively on the vortensity/entropy gradients within the planet's corotation region.For low-mass planets, M p 10M ⊕ , these two components of the gas torque are of the same order of magnitude but have opposite signs, leading to the so-called type-I migration regime (e.g., Tanaka et al. 2002).Positive corotation torques can dominate over the differential Lindblad torque locally -so that a planet migrates outwards transiently (e.g., Dittkrist et al. 2014;Baillié et al. 2016;Guilera et al. 2017).However, overall, planets tend to experience fast inward migration on a global scale.This makes it challenging to explain the existence of giant planets at moderate-to-large distances from the central star (e.g., Miguel et al. 2011;Ronco et al. 2017;Voelkel et al. 2022).If planetary cores of the order of 10M ⊕ can hover around these regions while there is still sufficient gas in the disk, this may lend a possible way to form giant planets at these locations.
In order to alleviate the fast-inward migration problem, different mechanisms have been proposed to slow down -or even revert -the migration of growing planets.Most of these focus on thermodynamic effects in non-isothermal disks or modifications of temperature/entropy and gas density gradients close to the planet (see Paardekooper et al. 2022, for a recent review, and references therein).In addition, Guilera & Sándor (2017) and Guilera et al. (2020) show that pressure maxima induced by a viscosity transition in the disk can act as a migration trap for growing planets.Moreover, Benítez-Llambay et al. (2015) show that the luminosity of an accreting planet is able to modify the temperature/density close to the planet asymmetrically in a way that an additional "heating torque" develops.The heating torque is always positive and can slow down or even revert the planet's migration.More recently, Guilera et al. (2019) and Guilera et al. (2021) incorporated this effect into a global model of planet formation, showing that heating torques can generate significant outward migration for planets growing by pebble accretion.
Given the fact that protoplanetary disks are made of gas and a small fraction of solids -elements condensed as grains or dust-it is relevant to consider the possibility that gravitational torques exerted by solids can affect planet migration.Perhaps because solids are thought to contribute with about only 1% of the total mass of the protoplanetary disk, there have only been a handful of studies about the exchange of angular momentum between the planets and the solid component of the disk (e.g., Capobianco et al. 2011;Chrenko et al. 2017;Benítez-Llambay & Pessah 2018;Chen & Lin 2018;Kanagawa 2019;Pierens et al. 2019;Regály 2020).In particular, the 2D hydrodynamical simulations from BLP18 show that for low-to intermediate-mass planets, an asymmetric dust-density distribution develops and exerts a net torque onto the planet.This "dust torque" is in general positive and its value depends on the Stokes number of solids, the dust-to-gas mass ratio, and the planet mass.It can be of the order (or much larger) than the gas torque and can locally slow down, halt, or revert planet migration for planetary cores with M p 10M ⊕ .It is then imperative to consider the cumulative effect of the dust torque when integrating the equations of motion of the embedded forming planets.
The aim of this work is to build on the torque measurements provided by BLP18 and investigate the global dynamics of low-mass planetary cores embedded in dusty disks in order to quantify the impact of the dust torque on their migration and formation history.The rest of the paper is organized as follows.In §2, we present the scope of our framework for computing the torque exerted on a planet embedded in a gaseous disk that can contain a dust component.Subsequently, we compute torque maps, characterizing the dependence of the total torque on a planet on the various physical parameters describing the planet-disk system and/or the dust component.These include the planet's location, the dustto-gas mass ratio, the stellar mass accretion rate, and the dust Stokes number distribution.We present the results for steady-state locally-isothermal disks (with a radial temperature gradient) and non-isothermal disks (with radial and vertical thermal structure) in §3 and §4, respectively.Based on these torque maps, in §5, we compute for the first time the formation tracks of lowmass planets embedded in dusty disks where the torque exerted by the dust component is considered.This allows us to quantify the role of dust in low-mass planetary migration considering a range of solid-accretion rates on single planets embedded in both steady-state and evolving disks.In §6 we outline the implications of our findings, together with the most important present caveats and possible directions to address them.Finally, in §7 we summarize the most important takeaways and conclusions from this first study.

FRAMEWORK FOR COMPUTING PLANETARY TORQUES
Computing the torque exerted on an embedded planet by the gas and dust component in a protoplanetary disk requires solving the associated multi-species hydrodynamic equations for a particular disk model, planetary mass, and a dust-size distribution.This is, in all generality, a daunting task.Here, we approach this problem in a simplified way by considering existing models for the gas torque responsible for type-I migration (Tanaka et al. 2002;Jiménez & Masset 2017) together with the measurements of the dust-to-gas torque ratio provided by BLP18 obtained for wide range of planetary masses and dust Stokes numbers for a single planet embedded in a standard α-disk model in steady-state.

Torque Calculation
We consider the total torque, Γ tot , acting on a planet as given by1 where Γ g and Γ d are the torques exerted by the gas and dust, respectively.For a given disk model, these torques are obtained as described below.

Gas Torque
For the locally-isothermal disks discussed in §3, we employ the gas torque model associated with type-I migration as derived by Tanaka et al. (2002) assuming a globally isothermal disk with a power-law density profile Here, Γ g is the sum of the differential Lindblad and corotation torques and Γ p stands for the reference torque evaluated at the radial location of the planet r p as given by with Σ g the gas surface density, q = M p /M the planetto-star mass ratio, and h = c s /(rΩ K ) the disk aspect ratio.Here, c s is the sound speed and Ω K = GM /r3 is the Keplerian frequency.The subscript "p" indicates quantities being evaluated at r = r p .
Using Eq. ( 2) to compute the gas torque for disks with radial temperature gradients is of course an approximation but this enables us to capture the characteristic dependencies of type-I migration with the underlying disk model in a simple way.Note that, under the stated assumptions, Γ p encapsulates all the radial dependencies of the physical parameters characterizing the gas torque in a dust-free disk.
For the non-isothermal disk models (i.e., with vertical thermal structure) presented in §4, we employ the gas torque models responsible for type-I migration derived by Jiménez & Masset (2017).The gas torque Γ g we consider is then again the sum of the differential Lindblad and corotation torques and is also proportional to Γ p defined in Eq. (3).In this case, the differential Lindblad torque depends mainly on the local gradients of the

Dust Torque
If the dust-to-gas mass ratio is small, we can neglect the dust back-reaction force onto the gas and Γ d scales linearly 3 with , i.e., We furthermore assume that, similarly to the models for the gas torque4 , the dust torque (for a given dust-to-gas mass ratio) is proportional to the reference torque, i.e., where Γ p and Γ 0 stand for the reference torque in Eq. ( 3) evaluated at the location of the planet r p and at the reference location r = r 0 , respectively.Combining Eqs.(4) Note-Values obtained for a reference radius r0 = 1 au and dust-to-gas mass ratio 0 = 0.01.
and ( 5), the dust torque can thus be written as The value of the ratio Γ d (r 0 , 0 )/Γ 0 is obtained from BLP18, where r 0 = 1 au and 0 = 0.01 (see their Fig. 2).The dust torque can then be finally expressed as where (Γ d /Γ 0 ) BLP18 is the dust torque measured from numerical simulations in BLP18.These values are listed in Table 1 for logarithmically spaced planetary masses and Stokes numbers, St, in the range 0.1 ≤ M p /M ⊕ ≤ 10 and 10 −2 ≤ St ≤ 1, respectively.
To obtain values of Γ d for planetary masses and Stokes numbers that are not tabulated, we use bilinear interpolation in logarithmic space.In Fig. 1, we show the interpolated dust torques from Table 1, which have been calculated considering r p = r 0 = 1 au and = 0 = 0.01.
In order to use Eq. ( 7) to model planets migrating and growing in disks with more complex structure, we implicitly assume that Γ d does not depend strongly on the gas dynamics or the state of motion of the planet.We note here that the dust torque is expected to be more sensitive to the disk dynamics for small Stokes numbers, for which the dust torque decreases.Also, the dust torque could be sensitive to the state of motion of the planet if the relative drift between solids and the planet is highly affected.This could happen, for example, if migration is very fast (with respect to the radial drift).
We also note that the accretion of dust, i.e., the removal of solids inside the planetary Hill sphere, can affect the magnitude of the dust torque (e.g., Regály 2020).The validity and robustness of these hypotheses need to be assessed through dedicated hydrodynamical studies.

LOCALLY ISOTHERMAL α-DISKS
We first consider a steady-state, locally isothermal αdisk model (Shakura & Sunyaev 1973) with a dust component characterized by a constant Stokes number throughout the disk as used by BLP18 and described in Weber et al. (2018).The α-parameter relates the accretion rate to the gaseous disk surface density via where Ṁ is the mass accretion rate onto the central star, ν = αc s h is the disk viscosity, and the aspect ratio, h, is considered to be constant.Under these assumptions, the background gas surface density and temperature are power laws in radius, with

Torque Maps for Fiducial, Isothermal Disk Model
In what follows, we will employ the term "torque map" to characterize the dependence of the total torque Γ tot on the various physical parameters describing the planet-disk system and/or the dust component.It is convenient to present normalized values Γ tot /Γ p using the reference torque Γ p , defined in Eq. (3).
In order to understand the scales involved in the problem at hand, and build intuition before considering more elaborate models, we first consider dusty disk models similar to those in BLP18.These correspond to steadystate, vertically-isothermal α-disks with constant aspect ratio h = 0.05 around a central star of 1 M with Ṁ = 10 −7 M /yr and α = 3 × 10 −3 .When the dust component is present, we set its constant gas-to-dust mass ratio to = 0.01.
The normalized torque map Γ tot /Γ p as a function of planet mass and Stokes number, computed at 1 au for the isothermal disk model considered, is shown in Fig. 2. The regions of the torque map corresponding to negative/positive values are associated with inward/outward migration.The two distinct regimes identified by BLP18 are evident; the gas-dominated regime for Stokes numbers St 0.05-0.2and planet masses M p 3 M ⊕ , and the gravity-dominated regime for St 0.3-0.5 and planet masses M p 7 M ⊕ .
In the gas-dominated regime (low values of St) the torque increases with the Stokes number, allowing the outward migration of more massive planets as St increases (up to a masses of 3 M ⊕ at St 0.1).On the other hand, in the gravity-dominated regime, the total torque is less sensitive to the value of St, and the maximum mass for outward migration stays at about M p 5 M ⊕ to 6 M ⊕ .In between these two regimes there exists a transition region (St 0.2-0.5)where the intensity of the dust torque is strongly reduced and can even change sign resulting mostly in inward migration.The total torque in this transition regime can depend sensitively on the mass of the planet.For the stan-dard value of the dust-to-gas mass ratio considered here, = 0.01, the total torque is practically equal to the gas torque for low Stokes numbers (St 0.02).However, as we will show in § 4.2 this situation can drastically change if the dust-to-gas mass ratio increases.
The normalized torque map Γ tot /Γ p as a function of planet mass and disk location associated with our fiducial disk model is shown in Fig. 3.The leftmost panel corresponds to a dust-free disk whereas the subsequent panels show the results for dusty disks with three different (constant) Stokes numbers St= {0.1, 0.3, 1}.For the dust-free disk, Γ tot = Γ g and Γ tot /Γ p is, by construction, constant with Γ tot /Γ p = −1.634.In the cases corresponding to St= 0.1 and St= 1, the dust torque generates a significant region of positive total torque, depicted with blue.For St= 0.1, this positive total torque region extends up to 3 M ⊕ , while for the case of St= 1 the positive total torque region extends up to 5 M ⊕ .The situation is quite different for the case of St= 0.3 where outward migration only ensues for masses 1 M ⊕ M p 6 M ⊕ and becomes inward outside that range.
Note that for the isothermal disk model under consideration, all radial dependencies are encapsulated in Γ p and the value of the total torque depends only on the mass on the planet and not on its location.Furthermore, because the disk aspect ratio is known, in this case h = 0.05, Γ p can be directly computed given the planet mass and planet location, and the total torque can be easily reconstructed from the color magnitude of the maps.In the next section, we will see that for non-isothermal disks (where the disk aspect ratio h is no longer constant), the transition between positive and negative total torques depends on the planet location; with implications for its migration history.
Having gained some insight into the relative contribution of the dust torque in the framework of the isothermal α-disk, we consider next a more realistic disk model.

NON-ISOTHERMAL α-DISKS
We consider a steady-state, non-isothermal axisymmetric α-disk with radial and vertical thermal structure.The disk model is obtained by solving the classical structure and transport equations in the vertical direction, considering that vertical hydrostatic equilibrium at each radius results from viscous heating and irradiation from the central star5 .Specifically, in order to compute the vertical structure of the disk, for a given α-viscosity parameter and stellar mass accretion rate Ṁ , we solve at each radial location the set of equations where P , ρ, Ω, T , F and z are the pressure, density, angular frequency, temperature, heat flux, and vertical distance from the disk midplane, respectively, and ∇ = d ln T /d ln P (see Guilera et al. 2017Guilera et al. , 2019, for further details).
The logarithmic gradients for the gas surface density and mid-plane temperature, as well as the disk aspect ratio h, for a non-isothermal disk with α = 3 × 10 −3 and Ṁ = 10 −7 M /yr are shown in Fig. 4. The (constant) values corresponding to the isothermal disk model are shown with dashed lines in the same figure for comparison.It is evident that the radial structure of the nonisothermal disk model cannot be described by a unique power law.In these models, heat transport and, in particular, radiative opacities play a key role in determining the density and mid-plane temperature profiles of the disk.The non-trivial radial disk structure may lead to noticeable differences between the non-isothermal and isothermal gas torque maps.However, in order to obtain the dust torque for the non-isothermal disk we assume that the dust torque can still be computed using Eq. ( 7), where the non-isothermal disk aspect ratio is used when computing the reference torque Γ p . .Radial profiles of the logarithmic gas surface density and mid-plane temperature gradients, and the disk aspect ratio, for the steady-state solution of an α-disk considering Ṁ = 10 −7 M /yr and α = 3×10 −3 .The solid lines correspond to the non-isothermal disk, while the dashed ones illustrate the isothermal disk case.The abrupt changes in the mid-plane temperature and gas surface density gradients are due to sharp disk opacity transitions at such locations.
The normalized torque map Γ tot /Γ p as a function of planet mass and Stokes number, computed at 1 au for the non-isothermal disk with a dust-to-gas mass ratio = 0.01, is shown in Fig. 5.As in the case of the isothermal disk model, the two regions of positive total torque, as well as the transition between these at St 0.2-0.5, are clearly visible.Similarly to the vertically isothermal disk, the dust torque does not play a relevant role for the standard dust-to-gas mass ratio of 0 = 0.01 when dust particles with low Stokes numbers e.g., St 0.02 are considered.All of these features in the torque maps evaluated at 1 au, as seen in Fig. 5, are also broadly present in these torque maps when calculated at other disk locations.However, due to the radial disk structure seen in Fig. 4, the value of Γ tot /Γ p does depend on the planet's location, in contrast to the isothermal disk case.

Dependence on Planet Location
The normalized torque map Γ tot /Γ p as a function of planet mass and location for the non-isothermal disk model with gas torque as provided in Jiménez & Masset (2017) is shown in the top row of Fig. 6.The leftmost panel corresponds to a dust-free disk whereas the subsequent panels show the results for dusty disks with three different (constant) Stokes numbers St= {0.1, 0.3, 1}.Similarly to the vertically isothermal disk, the total gas torque is always negative, implying inward planet migration.However, in this case the gas torque presents a slightly more complex dependence on the location of the planet (see below).
The results for the dusty disks displayed in Fig. 6 are qualitatively similar to the case of the isothermal disk: for the cases of St= 0.1 and St= 1, the dust torque generates a significant region of positive total torque.In these two cases, differently, than for the isothermal disk case, the value of the planet mass at which there is a transition between negative and positive torque regions depends on the planet's location in the disk, see Fig. 4.This non-trivial radial dependence observed also to a lesser extent in the dust-free disk is due to the fact that the radial structure for the non-isothermal disk deviates from power laws and that the Lindblad and corotation torques as provided in Jiménez & Masset (2017) depend on the disk diffusivity and viscosity.
For disks with dust particles with St= 0.3, the result is also qualitatively similar to the vertically isothermal disk.We find a total negative torque throughout the disk for low-mass (up to about 1 M ⊕ ) planets.As the mass of the planet increases (at a fixed location), the total torque becomes positive, but it becomes negative again beyond a certain mass.As explained above, for non-isothermal disks the mass threshold for the transition between regions of positive and negative total torque depends on the location of the planet due to the radial disk structure seen in Fig. 4.
In order to illustrate the relative difference between the torque maps obtained for isothermal and non-isothermal disk models considered, the bottom row of Fig. 6 shows the associated differential torque map.The aim of this figure is to show how the radial dependence between both models differs, especially when the dust torque is included in the computation of the total torque.The radial structure of the normalized differential torque map in the dust-free case simply reflects the structure of the normalized non-isothermal gas torque map (the leftmost panel on the top row in Fig. 6) because the normalized isothermal gas torque map is constant.When dust is present, the radial structure of the differential torque map is more complex.For St= 0.1 and St=1 there is a trend, for low-mass planets wherein the total torque is positive, the total isothermal torque is larger (lower) than the total non-isothermal torque within (beyond) ∼ 10 au.For larger planets wherein the total torque is negative the trend is the opposite, the total isothermal torque is lower (larger) than the total non-isothermal torque within (beyond) ∼ 10 au.For the case of St= 0.3 and the regions where the total torque is negative, the total isothermal torque is more (less) negative than the non-isothermal within (beyond) ∼ 10 au.On the other hand, in the region where the total torque is positive, the total isothermal torque is lower (larger) than the non-isothermal within (beyond) ∼ 10 au.

Dependence on Dust-to-gas Mass Ratio
The value of the dust-to-gas mass ratio can in principle depart significantly with respect to the fiducial value = 0.01 and it is thus necessary to understand how the total torque on a planet changes accordingly.The normalized torque maps Γ tot /Γ p for the non-isothermal disk as a function of planet mass and dust-to-gas mass ratio, , for four different Stokes numbers St= {0.01, 0.1, 0.3, 1}, are shown in Fig. 7.These torque maps correspond to a planet located at 1 au but the results obtained at other locations in the disk are qualitatively similar.
As seen in Fig. 7, planets with M p 0.5 M ⊕ embedded in a disk with dust of small Stokes number St 0.01 can be subject to strong positive torques for dust-to-  gas mass ratios 0.05.As the dust-to-gas mass ratio increases, the total torque becomes positive for planets with larger masses.Planets with masses in the range 0.5 M ⊕ to 5 M ⊕ can still experience a positive net torque provided the dust-to-gas mass ratio increases from 0.05 to 1.This dependence of the total torque on the dust-to-gas mass ratio may play an important role in the formation of planets by pebble accretion inside the water ice-line, where Stokes numbers tend to be small (Guilera et al. 2023, in preparation).
The morphology of the torque maps for St= 0.1 and St= 1 is qualitatively similar to the case with St = 0.01.However, particles with these Stokes numbers, respectively associated with the gas-and gravity-dominated regimes discussed by BLP18, can exert large positive torques at significantly lower dust-to-gas mass ratios.For instance, planets with M p 0.3 M ⊕ experience positive total torques for dust-to-gas mass ratios as low as 2 × 10 −3 .Planets of higher masses require higher dust densities to experience positive total torques but it is remarkable that planets as massive as 10 M ⊕ experience positive torques for dust-to-gas mass ratios as low as 2 × 10 −2 .This behavior may have important implications for the formation of planets by pebble accretion outside the water ice-line (see Guilera et al. 2023, in preparation).
For Stokes numbers in the transition region (St 0.3) the torque map is very different, with positive total torques arising only when 0.01 and M p 1 M ⊕ (at 1 au).Due to the strong dependence of the torques with the location of the planet in this regime (see Fig. 6) we expect the range of masses and dust-to-gas mass ratios for outward migration to be sensitive to the planet's distance to the central star.

Dependence on Mass Accretion Rate
Protoplanetary disks typically evolve on time-scales of the order of 1 to 10 Myr.During this evolution, the stellar mass accretion rate Ṁ drops from about 10 −7 M /yr for the youngest proto-stars to 10 −9 M /yr for the oldest ones.This reduction in the stellar mass accretion rate is associated with a lower disk gas surface density (Hartmann et al. 1998).This motivates the analysis of the dependence of the torque maps on the mass accretion rate Ṁ .
The normalized torque maps as a function of planet mass and location for a non-isothermal disk with Ṁ = 10 −9 M /yr for a dust-free disk and for a dusty disk with three different (constant) Stokes numbers St= {0.1, 0.3, 1} are shown in Fig. 8.Note that even though the normalized torques are similar in magnitude to the case with Ṁ = 10 −7 M /yr, the gas torque is about one order of magnitude lower for the lower accretion rate.This reflects the fact that the surface density of α-disk models in steady-state decreases with decreasing mass accretion rate Ṁ according to Eq. ( 8).
The normalized torque map for the dust-free disk is very similar to the case with Ṁ = 10 −7 M /yr and the total torque is always negative.For disks with dust particles with Stokes number St= 0.1 the dust torque remains dominant, leading to a positive net torque, for planets with masses up to roughly 1.5 M ⊕ and 4 M ⊕ .For Stokes number St= 1 the total torque remains positive for planets with higher masses, of the order of 3 M ⊕ to 6 M ⊕ depending on the planet location.We also note that in these last cases, despite the fact that they present similar trends than those observed for Ṁ = 10 −7 M /yr, the radial dependence of the transition between positive and negative torques is different.In contrast to the behavior seen in Fig. 6, for the higher accretion rate Ṁ = 10 −7 M /yr, here the dust torque dominates closer to the central star.As the disk evolves and its accretion rate and surface density decrease, dust particles with either St= 0.1 or St= 1 may be relevant for enabling low-mass planets in the inner disk regions to migrate outwards during the late stages of the disk evolution (see §6).
For St= 0.3, similarly to the case with the accretion rate Ṁ = 10 −7 M /yr, at a given disk location, there is a range of masses for which the total torque is negative.This mass range decreases from 1M ⊕ M p 6M ⊕ in the inner disk to 2M ⊕ M p 3M ⊕ in the outer disk.At lower accretion rates the total torque is negative for almost all planet masses beyond r p 7 au.

Dependence on the α-viscosity Parameter
Here we analyze the dependence of the torque maps on the α-viscosity parameter.We repeat the computations  in §4 considering the values α = 10 −2 and α = 10 −4 .We note that in order to ease the comparison, we adopt different accretion rates than in the case of α = 3 × 10 −3 .If we were to use the same value Ṁ = 10 −7 M /yr the disk gas surface density corresponding to α = 10 −4 would be much higher than for α = 3 × 10 −3 .This is due to the fact that when the viscosity is decreased the gas surface density needed to maintain a constant accretion rate throughout the disk must increase accordingly.The opposite occurs for α = 10 −2 .Thus, in order to achieve similar gas surface density profiles, we adopt Ṁ = 3 × 10 −7 M /yr for α = 10 −2 and Ṁ = 3 × 10 −9 M /yr for α = 10 −4 .We note that even though the gas surface density radial profiles are similar in these three cases, the corresponding mid-plane temperatures are higher for larger α values and the local mid-plane temperature gradients differ.
The normalized torque maps for α = 10 −4 and α = 10 −2 are shown in Fig. 9 and Fig. 10, respectively.The trends observed are similar to the case with α = 3 × 10 −3 .For the dust-free disk the total torque is always negative.For St= 0.1 and St=1, the dust toque becomes the dominant contribution of the total torque, generating significant regions of positive torque.In the transition regime, for St= 0.3, we obtain negative total torques for the less massive and more massive planets and positive total torques for intermediate masses.The transition between regimes presents a complex radial dependence, and in the case of α = 10 −4 the total torque is negative for all masses for planets located at r p 0.2 au.
It is worth noticing that the main features of the torque maps, including the magnitude of the normalized torques, are rather insensitive to the precise value of α, a highly unknown parameter in real protoplanetary disks.For α = 10 −2 , dust diffusion may have a significant effect, particularly for particles with small Stokes number.However, the contribution of such particles to the dust torque is less significant.Additional hydrodynamical calculations, including the effect of dust diffusion in highly turbulent disks, are required to characterize the dependence of dust torques on viscosity.In all previous sections, we have considered that the dust component can be characterized by a single Stokes number.Here, we consider dust with a distribution of Stokes numbers.For illustrative purposes, and to ease the comparison with the results in the previous sections, we adopt again a steady-state, non-isothermal α-disk with α = 3 × 10 −3 and Ṁ = 1 × 10 −7 M /yr and a dust-to-gas mass ratio = 0.01.

Impact of Dust-size Distribution
We briefly recall that for a mass distribution with dn/dm ∝ m −β , the total dust mass is ∝ m 2−β .Thus, depending on whether β < 2 or β > 2, most of the mass of the system is contained either in the larger or the smaller particles, respectively.In addition, it is easy to show that the size distribution is given by dn/da ∝ a 3β−2 .In the Epstein drag regime, the Stokes number is proportional to the particle size and thus dn/dSt ∝ St 3β−2 .
We consider three cases covering different characteristic possibilities: a distribution for which most of the mass is in small particles, with β = 5/2 and dn/dSt ∝ St −5.5 , a homogeneous distribution, with β = 2 and dn/dSt ∝ St, and a distribution for which most of the mass is in larger particles, with β = 3/2 and dn/dSt ∝ St −2.5 .For context, in the classical fragmentation equilibrium case β = 11/6 and thus dn/dSt ∝ dn/da ∝ a −3.5 (Dohnanyi 1969), implying particle distributions for which most of the mass is in larger particles.
We represent the distribution of Stokes numbers with 21 evenly spaced bins in logarithmic-scale between St= 0.01 and St= 1.The dust species i contributes with a dust-to-gas mass ratio i = ϕ i , such that i ϕ i = 1 and thus i i = .In Fig. 11 we show the histograms of the mass weights for the three distributions considered.Note that in the case of dn/dm ∝ m −5/2 the four smallest Stokes numbers in the distribution contribute with 75% of the mass, whereas for dn/dm ∝ m −3/2 the four largest Stokes numbers in the distribution contribute 75% of the mass.In the homogeneous distribution each species has the same ϕ i = 1/21.The total dust torque can be computed as the sum of the dust torque exerted by each individual dust species6 where Γ i d can be obtained from Eq. ( 7) by setting = i for each species.
The torque maps obtained for each distribution are shown in Fig. 12.In the case where most of the mass is in the particles with lower Stokes numbers, the torque map is similar to the gas torque (left panel in Fig. 6) and the solids do not play a significant role.On the other hand, for the homogeneous distribution and for the distribution where most of the mass is contained in large Stokes number particles (middle and lower panels in Fig. 12) the dust torque dominates and generates regions of positive torque for planets with masses up to 3 M ⊕ and 7 M ⊕ , respectively.These results suggest that dustgrowth models in which most of the mass remains in the form of large particles (e.g., Birnstiel et al. 2012) favor outward migration of low-mass planets provided these particles have St 0.1-1.As we show in Guilera et al. 2023, these results may have important implications for the role of the dust torque on the migration of low-mass planets growing by pebble accretion.

PLANET MIGRATION
Using the models for the gas and dust torques as described in §2.1, we can integrate the equations of motion for an embedded planet to obtain its migration history in the disk.In order to accomplish this, we implement a simple model in which the planet mass grows by accretion of pebbles with a fixed Stokes number.For simplicity and consistency, we evolve only one planet at a time.We consider initial conditions at r p = {1, 5, 10} au and the initial mass of the planet is 0.1 M ⊕ in all cases.The orbital integration is terminated when any of the following conditions is met: the planet grows to 10 M ⊕ , or the planet reaches the inner/outer disk edge located at 0.1/100 au.We note here that in general these conditions are fulfilled in less than ∼ 1 Myr for the cases of moderate and high pebble accretion rates (10 −5 M ⊕ /yr and 10 −4 M ⊕ /yr, respectively) and in less than 5 Myr for the case of a low pebble accretion rate (10 −6 M ⊕ /yr).

Steady-state Isothermal Disks Models
We first consider a steady-state, vertically isothermal α-disk model with α = 3 × 10 −3 and constant stellar accretion rate Ṁ = 10 −7 M /yr.When dust is present, the dust-to-gas mass ratio is = 0.01 and the Stokes number takes a fixed value with St= {0.1, 0.3, 1}.We consider three constant values for the accretion rate of solids onto the planet according to Ṁp = {10 −4 , 10 −5 , 10 −6 } M ⊕ /yr.The resulting formation tracks in dusty disks are shown in the left column in Fig. 13 with solid lines.For reference, the formation tracks in the dust-free disks, exhibiting inward migration, are shown with dashed lines.We remark here that, under our steady-state assumption, the high mass accretion rate assumed implies that the disk gas surface density remains high during the planet's formation and migration.In this regard, this choice provides a reference case for the early stages of a more realistic disk evolution model.We note that, in steady-state, the associated integrated disk mass between 0.1 au and 40 au is about 0.1 M , which corresponds to a gravitational stable disk under the classical Toomre's stability criterion.The total gaseous disk mass between 0.1 au and 100 au is ∼ 0.3 M , which is too massive for typical protoplanetary disks.However, when we consider disk models that evolve in time below ( § 5.3), we implicitly assume that the disk is initially massive and compact and that, as the mass accretion rate decreases, the disk mass decreases and it spreads viscously.
As seen in Fig. 13 For St= 0.3, planets migrate inwards at early times in all the cases, in agreement with the torque map showed in Fig. 3.For the case of a solid accretion rate of Ṁp = 10 −6 M ⊕ /yr, planets do not reach the needed mass ( 1 M ⊕ ) to revert the inward migration and they reach the inner disk at 0.1 au very fast with masses under 0.5 M ⊕ .For Ṁp = 10 −5 M ⊕ /yr, planets revert the inward migration and migrate outwards until they reach 6 M ⊕ .After reaching this mass their migration reverses again so that they migrate fast inwards reaching 10 M ⊕ at the inner disk edge.For the highest solid accretion rate we consider, Ṁp = 10 −4 M ⊕ /yr, the planets grow fast enough so that they do not migrate significantly before they reach 10 M ⊕ .Whether or not dust is present has a more significant impact on the formation tracks when the solid accretion rate is lower.Naturally, planets that reach M p 10 M ⊕ fast enough do not have the chance to spend significant time in the region of the torque maps where the dust contribution is important.The situation is different when the accretion of solids is slower.In particular, for Ṁp = 10 −5 or 10 −6 M ⊕ /yr, and dust with St= 0.1 and 1, planets can remain in disk regions where they are subject to positive torques for longer times, affecting their formation tracks with respect to the dust-free case significantly.

Steady-state Non-isothermal Disks Models
We now consider a steady-state, non-isothermal α-disk model with α = 3 × 10 −3 and constant stellar accretion rate Ṁ = 10 −7 M /yr and evolve an embedded planet sweeping the same set of initial conditions and physical parameters we used for the vertically-isothermal disk.The formation tracks in steady-state, non-isothermal dusty disks are shown in the middle column in Fig. 13 with solid lines.For reference, the formation tracks in the dust-free disks, exhibiting inward migration, are also shown with dashed lines.
Overall, the planet formation tracks for non-isothermal disks are qualitatively similar to the ones computed for vertically-isothermal disks.Some of the results worth highlighting are: i) In dusty non-isothermal disks, in contrast to the isothermal case, the transition from outward to inward (or inward to outward) planet migration does not occur at the same planet mass, in agreement with Fig. 6. ii) In dusty, non-isothermal disks with St= 0.1 none of the planets reach 100 au when the solid accretion rate is Ṁp = 10 −6 M ⊕ /yr.These planets reach 2 M ⊕ and reverse their initial outward migration somewhere between ∼ 50 au and ∼ 80 au, depending on their initial locations.They all reach the inner disk edge with roughly the same mass ∼ 5 M ⊕ .A similar behavior is observed for the case of St= 1 and Ṁp = 10 −5 M ⊕ /yr.iii) Similarly to the isothermal disks case, for high rates of accretion of solids, e.g., Ṁp = 10 −4 M ⊕ /yr, the planet formation tracks are rather insensitive to the presence of dust.This is especially true for the planets initially located at 5 au and 10 au, due to fast planet formation timescales.

Evolving Non-isothermal Disks Models
We consider here a simple model for the temporal evolution of a non-isothermal disk as the stellar accretion rate decreases in time.We adopt the same model considered in Bitsch et al. (2015)  This parameterizes a smooth transition between a mass accretion rate of 10 −7 M /yr and 10 −9 M /yr in 5 Myr.
As time advances and the mass accretion rate decreases, the gas surface density decreases accordingly.As in the steady-state disk models considered above, we set α = 3 × 10 −3 .We compute a dense grid of 50 steadystate α-disk models with equally-spaced accretion rates Ṁ between the initial and final accretion rates.We interpolate between these 50 models as needed to obtain a smooth disk evolution in time.For simplicity, we assume that the dust-to-gas mass ratio remains constant as the disk density evolves in time.This provides a conservative approach to assess the impact of dust in altering the planetary formation tracks with respect to the dust-free disks.
In the background of these time-dependent, nonisothermal disk models, we evolve an embedded planet sweeping the same set of initial conditions and physical parameters we used for the steady-state disk models.
The resulting formation tracks in dusty disks are shown in the right column in Fig. 13 with solid lines.For reference, the formation tracks in the dust-free disks, exhibiting inward migration, are shown with dashed lines.Note that for the highest pebble accretion rate considered, Ṁp = 10 −4 M ⊕ /yr, the formation tracks are similar to the ones in steady-state disks, irrespective of the Stokes number.This is due to the fact that formation timescales are shorter than migration timescales and planets form fast enough that the disk evolution does not affect the planet formation process.
For the cases with pebble accretion rates Ṁp = 10 −5 M ⊕ /yr, the planet formation tracks present some differences with respect to the steady-state, nonisothermal dust-free disk.However, these are not as noticeable as for the solid accretion rate Ṁp = 10 −6 M ⊕ /yr.This is particularly the case for Stokes numbers St= 0.1 and St= 1, for which planets in steadystate disks end at 0.1 au with 5 M ⊕ and at 100 au with 2 M ⊕ , respectively.In contrast, in evolving disks, all the planets reach 5 M ⊕ in 5 Myr in moderate-to-extended orbits, i.e., between 5 au and 30 au, depending on their initial locations.As the gas torque on the planet decreases, due to the decreasing gas surface density of the disk, the planets are allowed to remain in orbits with moderate to extended radii for longer times.Clearly, disk evolution plays an increasingly significant role the slower the planets grow.

Impact of the Dust-to-gas Mass Ratio
As stated in §2, in the absence of dust-feedback on the gas, the dust torque scales linearly with the dust-togas mass ratio.In order to explore the impact of the dust-to-gas mass ratio on the planet formation tracks we consider the set of values = {0.02,0.05, 0.1}.If protoplanetary disks have dust-to-gas mass ratio similar to the metallicity of their host protostars, values of up to = 0.03 could be found in disks around metalrich stars.Higher values of dust-to-gas mass ratios, in at least some disk regions, could be achieved via other mechanisms.For instance, the dust-to-gas mass ratio can increase in the inner disk regions due to dust drifting from larger radii (e.g., Drazkowska et al. 2016).Dustto-gas mass ratios near unity have been envisioned to be reached in preferential locations in protoplanetary disks (e.g., Drazkowska & Alibert 2017).As an example, we consider a fixed dust-to-gas mass ratio throughout the disk and compute the planet formation tracks for St= 0.1 assuming that the gaseous disk evolves in time according to Eq. ( 12).
The impact of the dust-to-gas mass ratio on a sample of evolutionary tracks is shown in Fig. 14.It is seen that for dust-to-gas mass ratios with 0.05 planets can migrate significantly outwards even for the highest solid accretion rate we considered, ṀP = 10 −4 M/yr.For the models with lower solid accretion rates, the effect of larger dust-to-gas mass ratios is increasingly dominant for outward planet migration, irrespective of Stokes number.As we show in a follow-up paper (Guilera et al. 2023, in preparation), enhanced dust-to-gas mass ratios in the inner disk regions, e.g., due to dust drifting from outer disk regions, can have a significant impact on the migration of growing planetary cores inside the water-ice line.

SUMMARY & DISCUSSION
While planet migration has been extensively studied in gaseous disks, the effect of solids on this process has not yet been addressed.Measurements of dust torques have been reported by BLP18, but the effect of these torques on the migration of planetary embryos has remained an open question.Aiming at filling this gap, in this work we quantify for the very first time the impact of the torque arising from solids on planet migration and assess its role on the formation tracks of planetary embryos embedded in classical protoplanetary disk models.
Based on the torque measurements of BLP18, we computed what we refer to as "torque maps" in order to quantify the importance of the dust torque (with respect to the gas torque) in terms of the disk thermodynamics, the planet location, the stellar mass-accretion rate, the dust-to-gas mass ratio, the level of turbulent viscosity, and different dust-size distributions.Subsequently, we investigated how the migration of growing planets is affected by disk thermodynamics, dust content, and global evolution.
In order to make progress in an inherently complex problem, we have made a number of simplifying assumptions and approximations.This approach enabled us to obtain the first quantitative estimates for the impact of the dust torque under a very wide range of physical conditions.Some of the effects we have found could become more or less prominent when models alleviating the caveats mentioned below are considered.Because of this, the theoretical formation path that a planetary core follows may depend on the details of the models employed.Nevertheless, it is worth stating some of the consequences associated with the broad trends that emerge from the results we obtained using standard assumptions to model the disk and its solid content.

Implications
Our findings have a series of implications for the processes involved in the formation and evolution of planetary cores in dusty disks.These, in turn, determine not only the formation timescales, final masses, and location of the planetary embryos but, potentially, also the final compositions of the resulting planets.
Low-mass planets migrate outwards beyond the water ice-line -Dust growth models suggest that most of the mass in solids is contained in particles in the high-end of the particle size distribution (e.g., Birnstiel et al. 2012;Stammler & Birnstiel 2022).In addition, Venturini et al. (2020a) and Drazkowska et al. (2021) showed that the mass-averaged Stokes number of the dust distribution is St ∼ 0.1 beyond the water ice-line (see their Fig. 1 and Fig. 4, respectively).Thus, the results presented in § 4.5 imply that forming planets of up to a few Earthmasses located beyond the water ice-line migrate outwards, unless the dust-to-gas mass ratio becomes lower than ∼ 10 −3 (see § 4.2).
Small-Stokes numbers could matter in dust-rich disks -Under a number of circumstances, the dust-to-gas mass ratio could increase significantly in some disk regions.This could occur in the inner part of the disk due to radial pebble drift from the outer disk (e.g., Drazkowska et al. 2016;Venturini et al. 2020b), near the water ice-line due to water vapor re-condensation (e.g., Drazkowska & Alibert 2017;Schneider & Bitsch 2021), or in a pressure bump where pebbles could accumulate (e.g., Pinilla et al. 2012;Zhu et al. 2012;Weber et al. 2018;Morbidelli 2020;Guilera et al. 2020;Chambers 2021;Jiang & Ormel 2023).This implies that if the dust-to-gas mass ratio reaches, for example, = 0.1, even particles with Stokes numbers as small as St= 0.01 can exert signif-icant torques if most of the mass is in those particles (Fig. 7).
Forming planets in extended orbits in dust-rich disks -Disks with large content of solids favor outward migration.This may offer a mechanism for the cores that start their formation in the inner disk region to end up located far from the central star, at tens of au (see Fig. 13 and Fig. 14).Moreover, these planets that migrate outwards, up to tens of au, could act as seeds for the formation of more massive planets in wide orbits.This could offer a natural path to explain the formation of these controversial objects, which are usually invoked as responsible for the ring and gap structures observed in protoplanetary disks by ALMA (e.g., Andrews 2020).

Current Caveats and Future Lines of Work
Here, we state and assess our most important hypothesis and mention future lines of work aimed at building more realistic models.More accurate and consistent models for the gas (and dust) torque could be obtained with dedicated, systematic numerical simulations considering more detailed disk models.
Precomputed prescriptions for type-I migration -Even though the disk models we consider have a non-trivial radial and thermodynamic structure, we adopted the torque prescriptions of Tanaka et al. (2002), derived for globally isothermal disks with power-law radial structure, and Jiménez & Masset (2017), derived for adiabatic disks.These prescriptions for the gas torque neglect migration feedback (see, e.g., Paardekooper 2014), an effect that is more pronounced in low-viscosity disks (see the review of Paardekooper et al. 2022).
Caveats inherited from BLP18 -By construction, our results contain all the caveats involved in computing the dust torques in BLP18.In particular, the numerical simulations involved are two-dimensional and thus insensitive to the vertical disk structure.More realistically, in three dimensions, the results are expected to depend on the vertical distribution of the solid component, which is sensitive to the level of vertical diffusion of solids, as regulated by the level of turbulence in the disk.Additional studies of dust torques properly considering the turbulent diffusion of solids are needed to better quantify the magnitude of dust torques.Moreover, simulations in BLP18 only considered disks with radial (−1/2) power-law gas-surface density structures.How the magnitude of dust torques is affected by the underlying gas-disk model is still uncertain, and could potentially affect some detailed aspects of our results.In particular, we assumed that dust torque does not depend on the gas torque (but see Appendix § A), which implicitly assumes that the underlying gas-disk model is not significantly important.This, however, may not be necessarily the case; more specifically, for solids that are highly coupled to the gas.Also, in BLP18 it is assumed that solids behave as a pressure-less fluid, which is reasonable for solids that are very well coupled to the gas.However, when solids are barely coupled, this approximation could break down.Furthermore, the dust torque measurements of BLP18 do not include the effect of planetary motion.Migration modifies the relative drift between solids and the planet and is expected to modify the shape and magnitude of the dustdistribution asymmetry developed close to the planet, which could either enhance or diminish dust torques.Moreover, if the planet migrates, it may develop eccentricity.Additional hydrodynamical simulations are needed in order to quantify the impact of eccentricity on dust torques.It is also worth recalling that the effective smoothing length used in the two-dimensional simulations in BLP2018 can affect the measured values of dust torques, as confirmed by Regály (2020).The dust torque in the gravity-dominated regime is robust because it arises from a large-scale asymmetry.However, in the gas-dominated regime the dust torque occurs due to a subtle asymmetry in the dust density distribution close to the planet (BLP2018).In order to minimize spurious measurements due to unresolved dynamics, BLP18 conservatively cut off the inner half of the Hill sphere to compute the torques.Nevertheless, three-dimensional simulations, where the smoothing length is not a free parameter, are needed to quantify more accurately the dust torque in different regimes.Thermodynamic effects could be important for dust torques.For example, the equation of state used could affect the gas dynamics and may have a non-negligible effect on the solid dynamics (see, e.g., Miranda & Rafikov 2019).Also, the effect of pebble accretion can change the dust asymmetry around the planet and hence the magnitude of the dust torque (see, e.g., Regály 2020), and generate heat that would lead to heating torques (e.g., Benítez-Llambay et al. 2015;Masset 2017), an effect that, for simplicity, we do not consider in this first work.Finally, the dust back-reaction force (feedback) is neglected in BLP18.This implicitly assumes that the density contrast of solids remains low throughout the disk and close to the planet.However, large density contrast of solids is observed in the numerical simulations in BLP18, which may have an impact on the dynamics of gas and solids close to the planet.Further work including the dust feedback is needed to assess the impact of this approximation.

Extrapolating BLP18's results to different disk models -
In order to extrapolate the results of BLP18 to other disk models we assumed a particular scaling of the dust torque with the disk parameters, see Eq. ( 7).However, the dust torque could scale additionally with other physical disk parameters that we do not consider here.In order to assess the robustness of our results, in Appendix A we assume an alternative model where the dust torque is proportional to the gas torque, see Eq. (A2), and compare the results.We show that our main results and conclusions are insensitive to the specific choice.
Pebble accretion, Stokes number/dust-size distributions -In order to compute the migration of the growing planetary embryos, we adopted constant pebble accretion rates and Stokes numbers (see § 5).A number of models usually adopt solid particles with constant Stokes number -or constant size-throughout the disk with pebble fluxes that decay in time, decreasing the pebble accretion rate on the planet (e.g., Lambrechts et al. 2019;Ogihara & Hori 2020).In more realistic planet formation models that include dust growth and evolution, both the pebble Stokes numbers and the pebble flux evolve in time, and throughout the disk, in a more complex way (e.g., Venturini et al. 2020a;Drazkowska et al. 2021).
Isolated planets -Our study assumes that a planet evolves in the disk in isolation.In order to extrapolate our results to disks containing multiple planets it would need to be assumed that migrating planets can be modeled independently of each other.This assumption will not hold in compact systems where the gravitational interaction between cores and planets could be significant.Also, multiple cores migrating are expected to exert strong perturbations on the dust surface density that can propagate from the outermost planet inwards via radial drift (see, e.g., Fig. 4 in BLP2018).The conditions for significant interference between multiple planets could in principle be assessed with targeted numerical simulations.In addition, when multiple planets grow simultaneously if one of them becomes massive enough to open a gap, or a partial gap, on the disk, the flux of pebbles to the planets in inner orbits can be halted (e.g., Lambrechts et al. 2014;Weber et al. 2018).This situation could also affect the magnitude of the dust torque for the inner planets.
Disk gaps induced by photoevaporation -The classical disk models we adopted neglected, for simplicity, the effects of photoevaporation due to the central star.In low-mass stars, X-rays are the dominant source of photoevaporation (see, e.g., Kunitomo et al. 2021).Venturini et al. (2020b) showed that X-ray photoevaporation can open a gap in the gas disk at a few au on a typical timescale of 1-2 Myr.This gap prevents the pebble flux from reaching the inner disk, significantly affecting pebble accretion onto the planets there.Thus, more realistic gas disk evolution models should also affect the computation of the dust torque.In a follow-up paper (Guilera et al. 2023, in preparation), we incorporate the computation of the dust torque in PLANETALP, our global framework for modeling planet formation, in order to study its effects on planet migration in more realistic settings for modeling planet formation.
Disk viscosity -Our work is based on the assumption that the disk remains laminar so we can describe the gaseous torques by simple prescriptions.Also, we used measurements of dust torques in BLP18 in laminar viscous disks.However, planets embedded in low-viscosity disks can induce vortices, which affect not only the gas (see e.g.Paardekooper 2014;McNally et al. 2019), but also the dust torque.Thus, caution is needed when extrapolating our results to low-viscosity or wind-driven disks.

TAKEAWAYS AND CONCLUSION
Here, we summarize the most relevant takeaways from our study emphasizing the results that we expect not to depend too sensitively on our main assumptions and approximations.

Main Takeaways
-The dust torque can contribute significantly to the total torque acting on a planetary embryo, M p 10M ⊕ , for standard values of the dust-to-gas mass ratio, i.e., = 0.01 (see Figs. 2, 5, and 6 ).
-The dust torque is positive in the vast majority of the parameter space spanned by 0.1 ≤ M p /M ⊕ ≤ 10 and 0.01 ≤ St ≤ 1 for = 0.01.The relative magnitude between the gas and dust torque components depends on the Stokes numbers of dust particles.
-There are two regimes where the dust torque dominates significantly for M p 10M ⊕ ; a "gas-dominated", St 0.1, and a "gravity-dominated" regime, St 0.5.These regions are associated, locally, with outward migration.This outcome is robust for a wide range of disk parameters such as the stellar mass accretion rate, the α-viscosity parameter, and/or the disk thermodynamics adopted.
-In the transition region between the previous two, i.e., St 0.2-0.5 and M p /M ⊕ ≤ 10, the value of the dust torque is sensitive to all the parameters involved (e.g., viscosity, stellar-accretion rate, the planet location) and the direction of migration depends on details.
-If there are disk regions with enhanced dust-to-gas mass ratios of order = 0.1, even particles with Stokes numbers as low as St 0.01 can exert a significant torque (see Fig. 7).At these higher dust concentrations, the torque exerted by larger dust particles can influence the dynamics of planets with masses larger than M p 10M ⊕ .
-The results outlined above do not depend sensitively on the stellar mass accretion rate or the α-viscosity parameter ( § 4.3 and § 4.4, respectively) or whether the dust particles are assumed to have a unique Stokes number or are distributed in a range.
-The properties of the torques as described above dictate the local dynamics of a planetary embryo in a dusty disk.Given that the relative magnitude of the gas and dust torque components evolves as the embryo migrates and grows, its long-term dynamics depends on the accretion rate onto the planet, the global disk structure and/or its long-term evolution.
-Even though the details of the evolutionary tracks are sensitive to various physical parameters, we find that under the set of assumptions usually invoked dust torques can have a significant impact, during significant fractions of the formation history of planetary embryos, by halting or reversing inward migration.The impact of the dust component naturally decreases for planetary embryos that grow too fast, e.g., Ṁp 10 −4 M ⊕ /yr, or disks that run out of material too quickly (see Fig. 13) -The most important implications of our findings can be summarized as follows: i) Low-mass planets migrate outwards beyond the water ice-line if most of the mass in solids is in larger particles there.ii) The torque due to solids with small Stokes numbers, St 0.01, can play a dominant role in disks with a moderately high dustto-gas mass ratio, 0.02 − 0.05 if most of the mass is in those particles.iii) Under a wide range of conditions, dust torques could enable low-mass planetary cores formed in the inner disk to migrate outwards and act as the seed for massive planets at distances of tens of au, which could be responsible for the large-scale structures observed by ALMA.

Conclusion
Our most important conclusion is that, under a wide range of conditions that are usually assumed to hold when modeling planet formation processes, dust torques can have an important effect on the evolution of planetary embryos.This may well be an important fundamental piece missing in the puzzle of planet migration and formation (Morbidelli & Raymond 2016).The detailed outcome when computing evolutionary tracks depends on the dust content and the global evolution of the disk.However, it is clear that the fastinward migration of low-mass planetary embryos, known to affect gaseous disks, can be significantly alleviated as a natural consequence of the torques exerted by dust.As a result, the long-term orbital evolution of planetary embryos can differ substantially from those obtained in dust-free disks.Given the importance of the role of dust and planet migration in planet formation, we hope that our findings will motivate further studies of the processes involved to develop a more complete and coherent framework for making progress under more astrophysically realistic assumptions.
We are grateful to the referee for a thoughtful report and useful suggestions that helped us improve the manuscript.OMG is partially supported by PICT 2018-0934 from ANPCyT, Argentina.OMG and M3B are partially supported by PIP-2971 from CONICET (Argentina) and by PICT 2020-03316 from Agencia I+D+i (Argentina).OMG acknowledges support by ANID, -Millennium Science Initiative Program -NCN19 171.OMG and M3B also thank Juan Ignacio Rodriguez from IALP for the computation managing resources of the Grupo de Astrofísica Planetaria de La Plata.P. B. L. acknowledges support from FONDECYT project 1231205.MEP gratefully acknowledges support from the Independent Research Fund Denmark via grant no.DFF 8021-00400B and the Institute for Advanced Study, where part of this work was done.

APPENDIX
A. AN ALTERNATIVE APPROACH FOR COMPUTING THE DUST TORQUE An alternative way of computing the dust torque from the hydrodynamical simulations of BL18 consists of assuming that the dust torque Γ d is proportional to the gas torque Γ g instead of the reference torque Γ p , as it is done in § 2. where Γ g (r 0 , 0 ) is the gas torque from BLP18 computed adopting 0 = 0.01 and r 0 = 1 au7 .The ratio (Γ d /Γ g ) BLP18 , as a function of planet mass, can be obtained8 from the information provided in Tab. 1. Finally, Γ g (r p , ) = Γ g (r p ) is computed for the corresponding gas disk model as described in § 2.1.1.
In order to understand the impact of the model adopted for the dust torque, Fig. 15 shows the normalized total torque for the case of a non-isothermal disk (see § 4) when the dust torque is assumed to be proportional to the gas torque instead of the reference torque Γ p , as it was previously done using Eq. ( 7).As in Fig. 6, the left panel in Fig. 15 shows the total torque over a planet in a dust-free disk, i.e., Γ tot = Γ g .The following three panels show the total torque Γ tot = Γ g + Γ d for the Stokes numbers St= {0.1, 0.3, 1}, considered to be constant throughout the disk, and assuming a gas-to-dust mass ratio = 0.01.
The results shown in Fig. 15 and Fig. 6 are qualitatively similar.The main difference is that in the latter the transitions between negative and positive torques in dusty disks occur at constant planet mass.This is because when the dust torque Γ d is proportional to the gas torque Γ g they both have the same radial dependence.The fact that this boundary is independent of the planet mass can be seen by noting that the zero torque condition Γ tot = 0 corresponds to 0 = (1 + Φ(M p , St, )) × Γ g , which is satisfied by the same value of the planetary mass M p regardless of the location of the planet, r p .
We also note that in Fig. 15 the total torque is lower compared to the total torque in Fig. 6, where the color scale is strongly saturated (note that Γ p is the same in both cases).This is also appreciable in Fig. 16, which shows the normalized total torque at 1 au (considering that the dust torque Γ d is given by Eq.A2) in a Stokes number -planet mass diagram as in Fig. 5.We remark here that now the gas-and gravity-dominated regimes, where the dust torque generates a total positive torque, appear as disconnected regions in this diagram.
In order to assess the impact of the two approaches we outlined to compute the dust torque, we plot in Fig. 17 the radial profile of the dust torque as obtained for three different planet masses, M p = {0.5, 1, 5}M ⊕ and for three different Stokes numbers St= {0.1, 0.3, 1} from left to right, respectively.In spite of the fact that the magnitude of the dust  torque (shown with solid lines) is in general slightly larger when using Eq. ( 7), its radial dependence is qualitative, and even roughly quantitative, similar independently of the approach used to compute it.

Figure 1 .
Figure 1.Normalized reference dust torque, Γ d /Γ0, for a range of planet masses and Stokes numbers for a verticallyisothermal disk as obtained from the simulations in BLP18, with rp = r0 = 1 au and = 0 = 0.01.The values shown correspond to the bilinear interpolation of data in Table 1. .gas surface density and mid-plane temperature while the corotation torque has four contributions, three of them associated with the disk radial gradients of vortensity, entropy, and temperature, and a fourth one accounting for the viscous production of vortensity, see Jiménez & Masset (2017) and Eqs.(26)-(28) and Eqs.(A.10)-(A.23) in Guilera et al. (2019) 2 .

Figure 2 .
Figure 2. Normalized torque map, Γtot/Γp, for a range of planet masses and Stokes numbers at 1 au for a verticallyisothermal disk.The color scale accentuates the regions where the total torque changes sign.Red/blue tones indicate negative/positive total torque, which locally implies inward/outward planet migration.

Figure 3 .
Figure 3. Normalized torque maps, Γtot/Γp, for a range of planet masses and locations in an isothermal disk.The left panel corresponds to a dust-free disk, whereas the other panels illustrate the results for a disk with dust-to-gas mass ratio = 0.01 and particles with constant Stokes number St= {0.1, 0.3, 1}.Since for isothermal α-disks all the radial dependencies of the torque are encapsulated in Γp, the normalized torques do not depend on the planet's location.
Figure 4. Radial profiles of the logarithmic gas surface density and mid-plane temperature gradients, and the disk aspect ratio, for the steady-state solution of an α-disk considering Ṁ = 10 −7 M /yr and α = 3×10 −3 .The solid lines correspond to the non-isothermal disk, while the dashed ones illustrate the isothermal disk case.The abrupt changes in the mid-plane temperature and gas surface density gradients are due to sharp disk opacity transitions at such locations.

Figure 5 .
Figure 5. Same as Fig. 2, but for the case of the nonishotermal disk.

Figure 6 .
Figure 6.Top row: normalized torque maps, Γtot/Γp, as a function of planet mass and location for the non-isothermal disk model with α = 3 × 10 −3 , Ṁ = 10 −7 M /yr, and dust-to-gas mass ratio = 0.01.The disk model is described in detail in §4.Bottom row: normalized differential torque maps between the isothermal and non-isothermal cases.

Figure 11 .
Figure 11.Histograms of the mass-weights ϕi corresponding to the species i for each of the three mass distributions considered in § 4.5.

Figure 12 .
Figure 12.Same as Fig. 6 but considering a Stokes number distribution between St= 0.01 and St= 1.The top panel shows the case where particles with the smaller size/Stokes number contribute the most to the mass of solids in the disk.The middle panel represents the case of a homogeneous mass (or Stokes number) distribution, while the bottom panel corresponds to the case where the mass of solids in the disk is dominated by the larger particle-size/Stokes numbers.

Figure 13 .
Figure 13.Formation tracks for planets initially located at rp = {1, 5, 10} au evolving in different disk models.The results corresponding to steady-state isothermal disks (see § 3), steady-state non-isothermal disks (see §4), and time-dependent nonisothermal disks (see §5.3) are shown in the left, center, and right set of panels.Each row of panels shows a given Stokes number in St= {0.1, 0.3, 1}.Solid lines represent tracks of planets in dusty disks, with a constant dust-to-gas mass ratio = 0.01, growing via pebble accretion with Ṁp = {10 −4 , 10 −5 , 10 −6 } M⊕/yr (violet, green, and cyan lines respectively).Dashed lines show the corresponding formation tracks in dust-free disks.

Figure 17 .
Figure 17.Radial profiles of the dust torque for different planet masses and different Stokes numbers.The solid lines represent the dust torque computed via Eq.(7), while the dashed lines correspond to the dust torque computed employing Eq. (A2).The black lines in the middle panel (the case of St= 0.3) highlight the fact that for a plant mass of 0.5 M⊕ the dust torques are negative in both approaches.