© The Authors 2025

1 Introduction

Understanding planet formation and evolution has remained a difficult endeavour as many physical processes act simultaneously on different length and timescales. The use of global models in planetary population syntheses helps assess the interaction between these processes (see e.g. the review by Burn & Mordasini 2024), but such global models rely on a good understanding of the underlying individual processes and on having simple prescriptions that still capture their physical nature. One fundamental ingredient of global models of planet formation and evolution is the evolution of the protoplanetary disc, as it impacts how planets grow and migrate. Past models usually assumed discs with alpha turbulent viscosities α_vis ≳ a few × 10^–3 (Bitsch et al. 2015; Liu et al. 2019; Emsenhuber et al. 2021; Drążkowska et al. 2023). However, observational constraints on turbulent line broadening from the gas emission of protoplanetary discs –or on vertical or radial diffusion of dust from the modelling of discs’ continuum emission– indicate that α_vis ≲ 10^–3 in the outer part of discs (at a few tens of astronomical units, see review by Pinte et al. 2023). These values of α_vis are not consistent with what is needed to explain observed stellar accretion rates by angular-momentum transport through turbulence (Hartmann et al. 1998; Emsenhuber et al. 2023). This and the lack of physical processes to sustain such high levels of turbulence throughout the disc are part of the reason why magnetohydrody-namic (MHD) winds have gained popularity as a way to transport angular momentum and drive disc evolution (Bai & Stone 2013; Bai 2013). Although a protoplanetary disc that is threaded by a magnetic field is susceptible to the magneto-rotational instability (MRI), detailed modelling including non-ideal MHD effects showed that MRI is suppressed in large parts of the disc from one to several tens of astronomical units (au). Turbulence in this region is thus mostly set by hydrodynamical instabilities and is expected to imply α_vis ≲ 10^–4 (see review by Lesur et al. 2023). This coincides with the region where planet formation is thought to take place.

With this in mind, it has become highly pertinent to investigate planet formation and evolution in MHD wind-driven discs at low effective turbulence α_vis (e.g. Ogihara et al. 2015, 2018; Kimmig et al. 2020; Kimura & Ikoma 2022; Wafflard-Fernandez & Lesur 2023; Lau et al. 2024). However, moving towards low turbulent viscosities comes with many challenges for planet formation and evolution models as the absence of viscous heating leads to the loss of convergence zones for migration (Emsenhuber et al. 2021) and the absence of diffusion leads to a saturation of the corotation torque, ultimately leading to faster, inward-only type I migration. The presence of an MHD wind promotes flatter or even inverted surface density profiles in the inner disc, which slows down or even reverses inward migration (Ogihara et al. 2015). Moving further towards low turbulent viscosities also promotes the dynamical corotation torque, which can reduce (Paardekooper 2014; Ogihara et al. 2017; Ndugu et al. 2021) or, conversely, also speed up type I migration (Paardekooper 2014; Pierens 2015).

In this work, we studied the importance of the dynamical corotation torque on type I migration in low-viscosity discs. In Section 2, we propose a simple description of the dynamical corotation torque for use in 1D disc models. In Section 3, we discuss our 2D hydrodynamic simulations of disc-planet interactions to investigate the effect of vortensity mixing in the planet’s horseshoe region due to viscosity, which provides the parame-terisations needed for our 1D model. In Section 4, we verify our model by direct comparison between 2D hydrodynamical simulations and 1D simulations. This enables us to assess the importance of the dynamical corotation torque on planet formation in Section 5. A discussion follows in Section 6, which notably includes an alternative parametrisation of the dynamical corotation torque for 1D models, followed by a conclusion.

2 Dynamical corotation torque: Preliminaries and memory timescale

Classically, planet migration is assessed by calculating the disc torque onto a planet held on a fixed orbit. In 2D non-magnetised disc models, the torque comprises¹ the Lindblad torque (Γ_lin, angular momentum exchange between the planet and its wakes in the disc) and the static corotation torque (Γ_c,static, angular momentum exchange between the planet and the disc in its horseshoe region; hereafter HS region). Several studies, including Paardekooper et al. (2011), have studied the static corotation torque for planets on fixed orbits. In particular, when the disc’s turbulent viscosity is high enough, the static corotation torque remains unsaturated, whereas it saturates at low viscosity, leaving only the Lindblad torque.

When the planet and the disc drift relatively to each other, the corotation torque actually features an additional component called dynamical corotation torque (Γ_c,dyn). The latter itself has two contributions (see review by Baruteau et al. 2016). One is due to gas just outside the planet’s horseshoe region, which executes a single U-turn with respect to the migrating planet. The torque due to this orbit-crossing flow scales with the planet’s migration rate and has a positive feedback on migration. The second contribution is from the gas that remains trapped in the librating region of the planet as the latter migrates. The torque due to this librating flow also scales with the planet’s migration rate, but it has a negative feedback on migration. Whether the dynamical corotation torque has a net positive or negative feedback on migration amounts to comparing the (inverse) vortensity of the librating and orbit-crossing flows (Masset & Papaloizou 2003). Figure 1 provides an illustration of the various flows around the planet and the associated torques.

Following Paardekooper (2014) and McNally et al. (2017), the dynamical corotation torque can be written as $${\text{Γ}}_{\text{c},\text{dyn}}=\frac{3}{4}{r}_{\text{P}}{\Omega }_{\text{P}}^3{\displaystyle {\int }_{y_s}^{x_{\text{s}}}\left(I_{v,\text{lib}}-I_{v,\text{cross}}\right)x^2\text{d}x,}$$(1)

where r_P and Ω_Ρ denote the planet’s orbital radius and its Keplerian frequency, and I_v,lib and I_v,cross give the inverse vortensity of the librating and orbit-crossing flows. The quantity I_v is defined as I_v = Σ/{(∇ × υ) · e_z}, with Σ the surface density of the disc, υ its velocity field, and e_z the unit vector along the vertical direction. The relative motion between the disc and the planet implies that the planet’s horseshoe region is asymmetric. This is reflected by the bounds in the above integral, which feature x_s, the radial half-width of the horseshoe region, and the quantity y_s = x_s + ṙ_Pτ_lib/2, with ṙ_P being the planet’s migration rate and τ_lib = 8πr_P/(3Ω_Ρx_s) the libration timescale (see Figure 2 in Paardekooper 2014). In the present work, we neglected the radial velocity of the background gas, in contrast, for instance, to McNally et al. (2017).

As recalled above, the dynamical corotation torque depends on the difference in vortensity between the librating and the orbit-crossing flows. In the following sections, we propose a simple descriptive model for estimating both vortensities, assuming an inwardly migrating low-mass planet. In particular, we propose a semi-analytical model for the vortensity of the librating flow.

Fig. 1 Sketch illustrating perturbed vortensity around planet’s vicinity (black circle) with relevant flows and torques highlighted: the inner and outer spiral density waves inducing the Lindblad torque, horseshoe U-turns ahead of and behind the planet leading to the static corotation torque, and the librating and orbit-crossing flows giving rise to the dynamical corotation torque.

2.1 Inverse vortensity of the orbit-crossing flow

For a planet that migrates inwards, the orbit-crossing flow originates at the radial distance r_P – x_s. For a low-mass planet of which the wakes shock the disc far enough from the planet not to influence the horseshoe region, which is the case for planet-to-star mass ratios of q_P ≪ h³(r_P) –where h is the disc’s aspect ratio (Lin & Papaloizou 1993)– the inverse vortensity of the orbit-crossing flow takes its local, unperturbed value. This results in $$I_{v,\text{cross}}=I_v\left(r_{\text{P}}-x_{\text{s}}\right)\approx I_v\left(r_{\text{P}}\right)\times \left[1-\frac{x_{\text{s}}}{r_{\text{P}}}\left(\frac{3}{2}-\beta \right)\right],$$(2)

where we assumed that the disc’s surface density Σ ∝ r^–β and that the disc’s velocity is purely Keplerian². The half-width of the planet’s horseshoe region is computed as $x_{\text{s}}=C{r}_{\text{P}}\sqrt{q_{\text{P}}/h\left(r_{\text{P}}\right)}$, with C being a constant factor that varies with the softening length of the planet’s gravitational potential (see Eq. (49) in Paardekooper et al. 2011).

2.2 Inverse vortensity of the librating flow

The gas vortensity satisfies an advection-diffusion equation with an additional source term in (∇Σ × ∇Ρ)/Σ³ (P denotes pressure) that represents baroclinic generation of vortensity (see e.g. Eq. (22) in Baruteau & Masset 2013). For a locally isothermal disc, the baroclinic source term is proportional to dT/dr × ∂Σ/∂φ and vortensity can be generated across the wakes of the planet (Casoli & Masset 2009; Ziampras et al. 2024b). In the absence of a baroclinic source term, which is the case for barotropic flows (such as in globally isothermal discs), vortensity evolves solely through advection-diffusion. The vortensity of the librating flow then depends on the migration history of the planet and how fast turbulent viscosity mixes the librating vortensity with that of the orbit-crossing flow. We introduce a memory timescale, τ_memory, which is a measure of the time it takes for the librating vortensity to adapt to the vortensity of the ambient gas. A planet that migrates in a disc with non-zero vortensity gradient carries the vortensity from its past locations. We assume that at a given time, t, the planet carries the vortensity it had at the time t – τ_memory, when it was located at a radial distance of r_P – ṙ_Pτ_memory. More specifically, we assume that $$I_{v,\text{lib}}=I_v\left(\min \left[r_{\text{P}}-{\stackrel{·}{r}}_{\text{P}}{\tau }_{\text{memory}},r_0\right]\right),$$(3)

where the minimum function min[] ensures that the vortensity cannot originate from further away than the planet’s initial location r₀. As mentioned already, in this work we limited ourselves to inward-migrating planets (ṙ_P < 0), and for outward migration (ṙ_P > 0) one should replace the minimum function in Eq. (3) with the maximum one. Introducing the migration timescale τ_mig such that ṙ_P = r_P /τ_mig, one obtains $$I_{v,\text{lib}}=I_v\left(\min \left[r_{\text{P}}\times \left(1-{\tau }_{\text{memory}}/{\tau }_{\text{mig}}\right),r_0\right]\right)$$(4) $$=I_v\left(r_{\text{P}}\times \min \left[1-\frac{{\tau }_{\text{memory}}}{{\tau }_{\text{mig}}},\frac{r_0}{r_{\text{P}}}\right]\right)$$(5) $$=I_v\left(r_{\text{P}}\right)\times {\left(\min \left[1-\frac{{\tau }_{\text{memory}}}{{\tau }_{\text{mig}}},\frac{r_0}{r_{\text{P}}}\right]\right)}^{3/2-\beta }.$$(6)

Viscous diffusion will take I_v,lib towards the local unperturbed vortensity I_v(r_P) on a timescale similar to the viscous timescale over the HS region ${\tau }_{\text{vis}}=x_{\text{s}}^2/v$. For high turbulent viscosity v, the viscous timescale is short compared to the migration timescale, and I_v,lib is close to the local unperturbed value I_v(r_P), making the dynamical corotation torque ineffective. For low turbulent viscosity, the viscous timescale across the HS region is long compared to the migration timescale, and I_v,lib deviates progressively from the local unperturbed value, so the dynamical corotation torque is expected to increase in magnitude. We note that in the limiting case of vanishing turbulent viscosity, the planet carries along the vortensity at its initial location. Further, it can be seen that for a disc with zero vortensity gradient (β = 3/2), the dynamical corotation torque vanishes. For a planet migrating up or down a vortensity gradient (migrating into regions with higher or lower vortensity), the dynamical corotation torque will slow down or speed up migration. We refer the reader to Paardekooper (2014) for an in-depth discussion of the different regimes of the dynamical corotation torque.

The mixing process of vortensity in the planet’s HS region is complex and still needs to be understood in detail. In our model, this process is modelled by a single parameter τ_memory. We expect τ_vis to pose an upper limit on τ_memory. However, this needs to be confirmed by measuring the vortensity evolution in 2D hydro-dynamical simulations. In the following, we present results of 2D hydrodynamical simulations in order to gain insight into the vortensity evolution around migrating planets.

3 2D hydrodynamical simulations

3.1 Numerical setup

We ran 2D hydrodynamical simulations of disc-planet interactions. This involves solving the continuity and Navier-Stokes equations, which are closed by an ideal, locally isothermal equation of state (see Masset (2002) for the detailed equations). We used FARGO3D (Benítez-Llambay & Masset 2016) along with the orbital advection algorithm of Masset (2000) to solve these equations. Calculations are performed on a cylindrical grid [r, ϕ] with r ∈ [0.4, 2.5]r₀ and ϕ ∈ [–π, π], where r₀, the planet’s initial orbital radius, defines the code’s unit of length. We used N_r = 782 cells linearly or logarithmically spaced along r, and Ν_ϕ = 2346 cells linearly spaced along ϕ. This choice of N_r allowed us to resolve the horseshoe region at the planet’s initial location by ten cells. For the boundary conditions at the grid’s radial edges, we applied a power-law extrapolation for the surface density, a Keplerian extrapolation for the azimuthal velocity, and an antisymmetric boundary condition for the radial velocity (thereby ensuring that the radial velocity cancels out at the grid’s radial edges). We also set wave-damping zones for r/r₀ ∈ [0.4, 0.52] ∪ [1.92, 2.5] to avoid reflections of the planet wakes (see de Val-Borro et al. 2006).

The gas surface density was initialised as a power-law, Σ₀(r) = Σ₀(r/r₀)^–β, with Σ₀ the initial surface density at the planet’s initial location r₀. We assume the disc evolution is driven by radial transport of angular momentum and is parameterised by the classical α-viscosity approach (Shakura & Sunyaev (1973); the alpha viscosity is taken uniform and denoted by α_vis). The disc’s aspect ratio is given by h(r) = h₀(r/r₀)^f with h₀ = h(r₀). In our fiducial simulations we adopted h₀ = 0.05 and f = 0.5. This implies a uniform background temperature, which avoids sharp vortensity perturbations close to the planet that would be produced by the baroclinic source term ∝ ∂_ϕΣ ∂_rT in the vortensity equation (see Section 2.2). Results of simulations with a non-uniform temperature profile are presented in Section 6.1, and additional globally isothermal simulations with h₀ = 0.07 are shown in Appendix A.

A freely moving planet was introduced at the beginning of our simulations. The planet-to-star mass ratio was chosen q_P = M_P/M_★ = 10^–5, which corresponds to M_P ≃ 3.3 M_⊕ (M_⊕ corresponds to the mass of the Earth) around a solar-mass star. This is about the mass of a super Earth, which is not expected to open a gap in our disc model as q_P/h³ ≪ 1. We used a smoothing length of 0.6H (r_P) for the gravitational potential of the planet (Müller et al. 2012). To avoid spurious inward acceleration of the inwardly migrating planet due to the fact that self-gravity is discarded, the force exerted by each grid cell on the planet involves the disc’s local surface density subtracted by its azimuthal average (Baruteau & Masset 2008; Benítez-Llambay & Masset 2016). This was done in our FARGO3D runs by activating the BM08 option in our code’s setup.

Since the star stays at the origin of the reference frame, indirect terms were included. These are the indirect term of the planet exerted on the disc and the indirect term of the disc exerted on the planet. Since the disc is not self-gravitating, the indirect term of the disc on itself was not included (Crida et al. 2022).

Our simulations aim to investigate how planet migration impacts the evolution of the vortensity of the librating and orbit-crossing flows for different values of α_vis, β, and Σ₀. The latter was varied according to the initial Toomre Q-parameter at the planet’s initial location: $$Q_0=\frac{h_0{r}_0{\Omega }^2\left(r_0\right)}{\pi G{\Sigma }_0}.$$(7)

Our default value for Q₀ is 8, for which we ran simulations with α_vis ∈ {10^–3, 10^–4, 10^–5} and $\beta \in \left\{\frac{1}{2},1,\frac{3}{2},2\right\}$. We note that in a locally isothermal disc, β < 3/2 implies that the dynamical corotation torque has a negative feedback on migration and therefore slows down inward migration, whereas for β > 3/2 the dynamical corotation torque has a positive feedback on migration and accelerates inward migration (Paardekooper 2014). Furthermore, we ran additional simulations with Q₀ ∈ {12, 16, 24, 32}, α_vis ∈ {10^–4, 10^–5}, and $\beta \in \left\{\frac{1}{2},1\right\}$ to gain more insight into the negative feedback effect of the dynamical corotation torque. Table 1 shows the full set of simulations presented in Section 3.2. Simulations were run for 3000 orbits or were stopped when the planet reached r_p ≲ 0.5r₀.

In the following, results of simulations are expressed in code units (c.u.), where the code’s unit of mass is the mass of the star, its unit of length is the initial orbital radius (r₀) of the planet, as already stated above, and its unit of time is Ω^–1(r₀).

3.2 Results

In the following sections, we discuss the results of the 2D simulations (left column of Figure 2). The 1D simulations are introduced and discussed later in Section 4.

3.2.1 Migration rates and orbital distances

First, we present the results of simulations with Q₀ = 8 for various values of α_vis and β in order to investigate the different migration behaviours predicted by theory (see Section 2). In the top two panels of the left column of Figure 2, we show the time evolution of the planet’s migration rate and orbital distance. For a disc with no vortensity gradient (β = 3/2) the migration rate slowly increases as the planet moves inward, which simply reflects the radial dependence of the Lindblad torque since both the static and dynamical corotation torques vanish. The planet reaches the inner edge of the grid in ≳2000 orbits.

For negative vortensity gradients (β = 1/2 and β = 1), migration rates take smaller values than for β = 3/2 since the static and dynamical corotation torques are both positive. We note that decreasing α_vis from 10^–3 to 10^–4 leads to slightly higher migration rates, and then going further down to α_vis = 10^–5 reduces migration rates more significantly. This behaviour reflects how the static and dynamical parts of the corotation torque vary with viscosity, since decreasing viscosity decreases the static part but increases the dynamical part. We remind the reader that for the range of (low) viscosities probed in our simulations, the Lind-blad torque is independent of viscosity (Muto & Inutsuka 2009). Only at low α_vis and low β does the migration rate increase with time, more appreciably at the beginning of the simulations, highlighting the negative feedback of the dynamical corotation torque on migration.

Finally, for a positive vortensity gradient (β = 2), we observe runaway migration at a rate that increases with decreasing viscosity. This behaviour cannot be attributed to the static corotation torque, although it is negative here, since its magnitude decreases with decreasing viscosity. It is therefore due to the dynamical corotation torque. In these simulations, the planet reaches the grid’s inner edge in 1000–1500 orbits. This all reflects the theoretical expectations.

3.2.2 Vortensity of the orbit-crossing and librating flows

We are interested in the time evolution of I_v,lib – I_v,cross as it is what sets the dynamical corotation torque (see Eq. (1)). To measure I_v,lib and I_v,cross in our 2D simulations, we used a very similar approach to that of Wafflard-Fernandez & Baruteau (2020). For the orbit-crossing flow, since in our models the horseshoe U-turn timescale is short in comparison to the viscous diffusion timescale over the HS region, I_v,cross can simply be measured at the orbital radius of the planet and just behind it in azimuth (as the planet migrates inwards). In practice, it was averaged over ϕ – ϕ_p ∈ [–0.2, –0.1] at r = r_p (ϕ_p denotes the planet azimuth, which is zero in our simulations since our runs were carried out in the frame corotating with the planet). For the librating flow, two strategies are employed. When β < 3/2 (no runaway), I_v,lib is averaged over azimuth at the planet’s orbital radius, excluding a small region around the planet; in practice, it was averaged over ϕ – ϕ_p ∈ [–π, –0.4] ∪ [0.4, π] at r = r_p. Finally, when the planet goes into runaway migration (β = 2), the librating region is heavily distorted, leaving just an island of librating material in front of the planet in the azimuthal direction (as the planet moves inward). In that case, again in practice I_v,lib was averaged over ϕ – ϕ_p ∈ [0.4, 0.6] at r = r_p.

The upper panels in Figure 3 illustrate the above method for a simulation with β = 1, in which case we see that the evolution of the vortensities is well captured throughout the simulation. Note the orbit crossing flow originating from ≃r_P – x_s which is in agreement with Eq. (2). In the lower panel of the figure, we stress that, although both vortensities increase with time, the ratio I_v,lib/I_v,cross reaches a quasi-steady state after ~1500 orbits. The ratio is not necessarily expected to remain at a constant value as the migration rate and the width of the HS region vary with orbital distance.

Back to Figure 2, the lower left panel shows the time evolution of I_v,lib/I_v,cross in the simulations. For β < 3/2, the planet moves up the disc’s vortensity gradient and thus has a vortensity deficit in its HS region, which explains why I_v,lib/I_v,cross > 1. Quite remarkably, in all our simulations with β < 3/2, I_v,lib/I_v,cross reaches a near-stationary value, while the planets keeps migrating, after a characteristic timescale that depends on β and α_vis. The lower the viscosity, the longer it takes to reach a steady state, the higher the final value for I_v,lib/I_v,cross and the stronger the negative feedback of the dynamical corotation torque on migration, as already shown above. For the disc and planet parameters of these simulations, the dynamical corotation torque really starts impacting migration for α_vis ≲ 10^–4. For β > 3/2, the planet moves down the vortensity gradient and thus builds up a vortensity excess in its HS region, implying that I_v,lib/I_v,cross < 1. This time, the runaway implies that I_v,lib/I_v,cross keeps decreasing with time and does not reach a steady state value. Finally, for β = 3 /2 we see that I_v,lib/I_v,cross stays at unity regardless of turbulent viscosity, as expected.

Fig. 2 From top to bottom: time evolution of planet’s migration rate, orbital distance, and inverse vortensity ratio Iv,lib/Iv,cross for simulations with Q0 = 8. The results of our 2D (with FARGO3D) and 1D (with the Bern model) simulations are shown side by side for easy direct comparison. Lines are coloured by αvis and different line styles indicate varying slopes of the surface density, β. The migration rates of the 2D simulations (top left panel) are shown as moving averages over 100 orbits. For the 1D simulations (right), we additionally show results of simulations that discard the dynamical corotation torque as transparent lines. For an easier direct comparison of the most relevant simulations, see the left column in Figure B.1 of Appendix B.

Fig. 3 Evolution of disc vortensity for simulation with Q0 = 8, β = 1 and αvis = 10–5. The top panels show snapshots of the vortensity field around the planet at 500, 1500, and 2500 orbits. The vortensity deficit in the HS region, and the vortensity advection of the orbit-crossing flow from the inner to the outer edges of the HS region are clearly visible. The vortensities of the orbit-crossing and librating flows measured by the method described in Section 3.2.2 are shown as contour lines in pink and cyan, respectively. The grey vertical lines mark the orbital radius of the planet and the extent of its horseshoe region (see Section 2.1 for the expression of the horseshoe half-width xs). The bottom panel shows the time evolution of the measured vortensities for the orbit-crossing and librating flows. Crosses mark the results shown in the top panels. The time evolution of Iv,lib/Iv,cross is shown in black.

3.3 Estimating the memory timescale

As shown in Section 3.2.2, I_v,lib/I_v,cross eventually reaches a quasi-steady state when migration does not run away (β < 3/2). While for simulations with α_vis ≳ 10^–4 the final value of I_v,lib/I_v,cross can essentially be attributed to the background vortensity gradient, for simulations at lower viscosity (α_vis = 10^–5) I_v,lib/I_v,cross takes much larger values, which points to a memory effect for the vortensity of the librating flow (see Section 2.2).

To estimate the memory timescale in our 2D simulations, we simply made use of Eq. (6), assuming that 1 – τ_memory/τ_mig < r₀/r_P for I_v,lib/I_v,cross to reach a constant value. Reversing Eq. (6) then gives $${\tau }_{\text{memory}}={\tau }_{\text{mig}}\times \left(1-{\left(\frac{I_{v,\text{lib}}}{I_v\left(r_{\text{P}}\right)}\right)}^{\frac{1}{3/2-\beta }}\right),$$(8)

where we recall that the migration timescale τ_mig = r_P / ṙ_P. Although migration rates are close to stationary in simulations with β < 3/2 (see Fig. 2), the migration timescale varies over time due to the planet moving inwards, and so does the memory timescale.

We computed the memory timescale in an extended set of 2D simulations that include higher values for Q₀ (12, 16, 24 and 32) for $\beta \in \left\{\frac{1}{2},1\right\}$, and for α_vis ∈ {10^–4, 10^–5}. Results are shown in Figure 4. The left panel first shows the time evolution of I_v,lib/I_v,cross for α_vis = 10^–4. Two trends emerge: (i) for a given Q₀, I_v,lib/I_v,cross is larger for smaller β; and (ii) for a given β, I_v,lib/I_v,cross decreases with increasing Q₀. The viscous and libration timescales over the HS region happen to have similar values within about 30% (τ_vis ≃ 130 orbital periods at r₀, τ_lib ≃ 100 orbital periods at r₀). This implies that turbulent viscosity takes the vortensity of the librating flow towards the local unperturbed value on a timescale similar to the libration timescale and thus no clear memory effect is expected to occur. We thus suspect that the ratio I_v,lib/I_v,cross is essentially explained by the background vortensity gradient (see Eq. (2) with I_v,lib ~ I_v(r_p)). Evaluating Eq. (2) at the planet’s initial location actually gives I_v(r_p)/I_v,cross ≈ 1.014 for β = 1/2 and ≈ 1.007 for β = 1, which is in decent agreement (albeit slightly above) with the values obtained in the left panel of Figure 4.

The middle panel of Figure 4 now shows I_v,lib/I_v,cross for α_vis = 10^–5. Here, the viscous timescale is one order of magnitude larger than the libration timescale, and we observe a clear build-up of vortensity memory that reaches a quasi-steady state for sufficiently fast migration (Q ≲ 24). Having previously shown that the background vortensity gradient also contributes to the final value for I_v,lib/I_v,cross, we subtracted the contribution of the background vortensity gradient when calculating the memory timescale using Eq. (8). The right panel of Figure 4 shows the obtained memory timescale divided by the local viscous timescale over the HS region (half-width). Although simulations with different β values reach a quasi-steady state for different values of I_v,lib/I_v,cross, they show very similar memory timescales. The panel also shows that the memory timescale has a mild dependence on the migration rate (via Q₀), but overall τ_memory ~ [0.3–0.4] × τ_νis is a good match to our simulations results. These experiments allow us to conclude that two effects impact the ratio I_v,lib/I_v,cross: the background vortensity gradient and a memory effect at low viscosity for the vortensity of the librating flow, which acts on a typical timescale ≃0.4τ_vis.

Fig. 4 Left and middle panels: time evolution of Iv,lib/Iv,cross in simulations with varying Q0 ∈ {8, 12, 16, 24, 32}, $\beta \in \left\{\frac{1}{2},1\right\}$, for αvis = 10–4 (left) and αvis = 10–5 (middle). Right panel: time evolution of ratio between memory timescale and viscous timescale across HS region for αvis = 10–5 simulations.

4 1D simulations

In this section, we present the results of 1D simulations that model the disc evolution and the planet’s orbital evolution. We compare them with the results of the 2D simulations of the previous section in order to test our new formulation for the dynamical corotation torque based on a memory timescale argument. The 1D simulations use parts of the Bern model of planet formation and evolution (Alibert et al. 2005; Mordasini et al. 2012; Emsenhuber et al. 2021).

4.1 Numerical setup

The 1D simulations model the evolution of a protoplanetary disc by solving the classical viscous diffusion equation for the disc’s surface density: $$\frac{\partial \Sigma }{\partial t}=\frac{1}{r}\frac{\partial }{\partial r}\left[\frac{3}{r\Omega }\frac{\partial }{\partial r}\left(r^2\Sigma {\alpha }_{\text{vis}}{c}_{\text{s}}^2\right)\right],$$(9)

with c_s = hrΩ being the sound speed and α_vis the alpha turbulent viscosity. As in our 2D simulations, the disc’s aspect ratio was chosen h = h₀(r/r₀)^f with f = 1/2 (uniform temperature profile), the initial surface density was Σ₀(r/r₀)^–β, and α_vis was taken constant. The 1D grid used 1244 cells logarithmically spaced over r ∈ [0.2, 2.5] r₀. For the boundary conditions, we used a power-law extrapolation again for the surface density. The default time step was one orbit at the planet’s initial location with additional constraints to ensure that the Courant-Friedrichs-Lewy condition is fulfilled and the planet does not move more than 10^–4r_P per time step.

At the beginning of the simulations, a planet was inserted at r₀ with the same planet-to-star mass ratio as in our 2D simulations. The planet’s semi-major axis was evolved according to the disc torques. For the Lindblad torque and the static corotation torque, we used the torque formulae of Paardekooper et al. (2011)³. For the dynamical corotation torque, Eq. (1) gives $${\Gamma }_{\text{c,dyn}}=\frac{1}{4}{r}_{\text{P}}{\Omega }_{\text{P}}^3\left(I_{v,\text{lib}}-I_{v,\text{cross}}\right)\left(x_{\text{s}}^3-y_s^3\right),$$(10)

assuming I_v,lib and I_v,cross are uniform over the HS region. We further set y_s such that $$y_s=\max \left(x_{\text{s}}+{\stackrel{·}{r}}_{\text{P}}{\tau }_{\text{lib}}/2,0\right),$$(11)

which ensures that in the fast migrating regime (runaway), the dynamical torque reaches $${\Gamma }_{\text{c,dyn}}=\frac{1}{4}{r}_{\text{P}}{\Omega }_{\text{P}}^3\left(I_{v,\text{lib}}-I_{v,\text{cross}}\right)x_{\text{s}}^3$$(12)

(as discussed in McNally et al. (2018); see their Section 6.2). This makes the calculation of the dynamical corotation torque valid in both the slow and fast migration regimes, thus independently of β. Finally, I_v,cross was calculated using Eq. (2) and I_v,lib via Eq. (6) with $${\tau }_{\text{memory}}=0.4{\tau }_{\text{vis}}\times \max \left(1-\frac{{\tau }_{\text{lib}}}{2{\tau }_{\text{vis}}},0\right).$$(13)

Eq. (13) ensures that a memory effect only sets in if the viscous timescale is sufficiently long compared to the libration timescale. If viscosity is indeed able to equilibrate the vortensity of gas that goes from the front to the back of the planet in the HS region, which takes about half a libration timescale, there will be no memory effect. A memory effect sets in at τ_vis ≃ 2τ_lib, corresponding to α_vis ≃ 10^–4 for our fiducial values (q = 10^–5 and h = 0.05). This coincides with the point where the static corotation torque saturates (see Fig. 2 in Paardekooper et al. 2011 and Fig. 11 in Baruteau & Masset 2013). The maximum function in the above equation ensures that τ_memory cannot be negative if the viscous timescale is shorter than half the librating timescale.

4.2 Results

Our results of 1D calculations are displayed in the right column of panels in Figure 2 for Q₀ = 8 and the same range of values for β (1/2, 1) and α_vis (10^–5, 10^–4, 10^–3) as in the 2D simulations. Overall, our results of 1D calculations are in very good agreement with those of the 2D simulations. In the calculations with no vortensity gradient (β = 3/2), the disc torque is entirely given by the Lindblad torque and does not depend on α_vis. We note that the 2D torque starts to slightly exceed the 1D torque over long timescales (≳1500 orbits), which is probably due to the gradual onset of non-linearities in the 2D simulations brought about by the planet’s wakes.

In discs with a negative vortensity gradient (β < 3/2), it is interesting to notice the initial stage where the migration rate decreases in the simulations with low viscosity (α_vis = 10^–5). It arises from the gradual increase in the dynamical corotation torque as a vortensity deficit builds up in the planet’s librating region. This stage is over once the planet has migrated far enough that the memory of the initial vortensity is lost (r_P – ṙ_Pτ_memory < r₀). A similar behaviour can actually be seen for the 2D simulations, which we can now attribute to the gradual increase in the I_v,lib/I_v,cross ratio. For higher turbulent viscosities (α_vis ≳ 10^–4), I_v,lib/I_v,cross is almost entirely set by the background vortensity gradient as seen in Section 3.3. However, Eq. (2) tends to somewhat overestimate I_v,lib/I_v,cross compared to the 2D simulations (compare the bottom panels), which leads to slightly lower migration rates, and therefore slightly slower migration in the 1D calculations (compare the top and middle panels). Note that we are able to reproduce with the 1D model the higher migration rates obtained for α_vis = 10^–4 compared to α_vis = 10^–3, discussed in Section 3.2.1. This indicates that the analytical timescale picture indeed captures the governing physical processes.

Now turning our attention to discs with a positive vortensity gradient (β > 3/2), we see that good agreement is obtained between the 1D and 2D simulations for α_vis ≥ 10^–4. For our lowest viscosity (α_vis = 10^–5), our 1D model tends to underestimate I_v,lib/I_v,cross, thereby predicting slightly faster runaway migration than in our 2D simulations. Yet, the agreement remains satisfactory and shows that we are able to model runaway type I migration driven by the dynamical corotation torque.

Before concluding this section, we point out that additional 1D calculations were carried out by discarding the dynamical corotation torque. The results of these 1D calculations, in which the disc torque only comprises the Lindblad and static corotation torques, are shown as transparent curves in the right panels of Figure 2. These extra simulations highlight the importance of including the dynamical corotation torque for the disc models we considered. More generally, our results show that the dynamical corotation torque must be taken into account in models of low-mass planet formation, which we further emphasise in the next section.

5 The extent to which the dynamical corotation torque reduces the type I migration torque

The two previous sections demonstrate that our new formulation for the dynamical corotation torque, which is based on a memory timescale for the evolution of the librating flow’s vortensity in the planet’s horseshoe region, reproduces well the migration behaviour of a low-mass planet obtained in 2D hydrodynamical simulations. Eqs. (6) and (13) show that the amplitude of the dynamical corotation torque depends on how the viscous timescale across the HS region (τ_vis) compares with the migration timescale (τ_mig) and the libration timescale (τ_lib) through the ratios τ_vis/τ_mig and τ_vis/τ_lib. For a given disc model (density, temperature, alpha viscosity), the amplitude of the dynamical corotation torque thus depends on the planet’s mass and orbital distance. To assess how the dynamical corotation torque impacts type I migration, we could carry out a number of 1D simulations (such as those of the previous section) with varying planet masses and initial locations, possibly by modelling the mass growth of the planet. For the present work, we adopted a simpler, more pragmatic approach by providing a simple estimate for how much the dynamical corotation torque is expected to reduce the classical type I migration torque, which is exclusively comprised of the Lindblad and static corotation torques.

We adopted a locally isothermal, low-viscosity disc model with shallow radial profiles of density and temperature, such that the dynamical corotation torque slows down type I migration. We shall assume the migration timescale is long compared to the libration timescale (i.e. $\stackrel{·}{r_{\text{P}}}\ll 3{\Omega }_{\text{P}}{x}_{\text{s}}^2/\left\{4\pi r_{\text{P}}\right\}$) so that Eq. (1) can be recast as (Paardekooper 2014) $${\Gamma }_{\text{c,dyn}}=2\pi \left(1-\frac{I_{v,\text{lib}}}{I_{v,\text{cross}}}\right){\Sigma }_{\text{P}}{r}_{\text{P}}^2{x}_{\text{s}}{\Omega }_{\text{P}}{\stackrel{·}{r}}_{\text{P}}.$$(14)

Writing the total disc torque Γ as Γ_lin + Γ_c,static + Γ_c,dyn and as $\frac{1}{2}$ M_Pr_PΩ_Pṙ_P (neglecting the time evolution of the planet’s mass), one arrives at $$\Gamma =\frac{{\Gamma }_{\text{lin}}+{\Gamma }_{\text{c,static}}}{1-\frac{M_{\text{HS}}}{M_{\text{P}}}\left(1-\frac{I_{v,\text{lib}}}{I_{v,\text{cross}}}\right)},$$(15)

where M_HS = 4πr_Px_sΣ_P is the mass inside the HS region. The total disc torque on the planet can thus be expressed as a simple correction to the classical type I torque, $$\Gamma =f_{\text{red}}\times \left({\Gamma }_{\text{lin}}+{\Gamma }_{\text{c,static}}\right),$$(16)

with f_red < 1 being a reduction factor that reads $$f_{\text{red}}={\left[1-\frac{M_{\text{HS}}}{M_{\text{P}}}\left(1-\frac{I_{v,\text{lib}}}{I_{v,\text{cross}}}\right)\right]}^{-1}.$$(17)

When Γ_c,dyn → 0, f_red → 1, whereas when Γ_c,dyn is large, f_red ≪ 1. Again, for a given disc model, f_red depends on the planet’s mass and orbital radius.

We evaluated the reduction factor f_red for disc models similar to those used in global models of planet formation and evolution: $$\Sigma \left(r\right)={\Sigma }_{5\text{au}}\times {\left(r/5\text{au}\right)}^{-\beta ,}$$(18) $$T\left(r\right)=T_{1\text{au}}\times {\left(r/1\text{au}\right)}^{-1/2},$$(19)

with Σ₅ au = 300 g cm^–2, β ∈ [1/2, 1], and T₁ au = 280 K. Our temperature profile (decreasing as r^–1/2) mimics that of a passively irradiated disc. To compute the dynamical corotation torque, we first computed the classical type I migration torque, which we dub the static torque, via the torque formulae in Paardekooper et al. (2011), and we used the migration rate resulting from the sole static torque to compute the dynamical corotation torque through Eq. (14). As in Section 4, I_v,cross is obtained with Eq. (2) and I_v,lib via Eqs. (6) and (13), assuming the planet has migrated sufficiently far for it to lose its memory of the initial location (i.e. r_P – r₀ < ṙ_P · τ_memory). We then iterated, using the migration rate associated with the total torque, to recompute the dynamical corotation torque until the value converged. The reduction factor is finally calculated as 1 + Γ_c,dyn/(Γ_lin + Γ_c,static). We note that, although our memory timescale formulation for the dynamical corotation has been worked out for uniform temperature profiles, it also reproduces the results of 2D hydrodynamical simulations with a non-uniform temperature profile well (locally isothermal discs), as we show in Section 6.1.

Figure 5 displays the reduction factor f_red as a function of planet mass and orbital distance for four disc models with α_vis ∈ [10^–5, 10^–4] and β ∈ [1/2, 1]. Each panel shows two regimes, detailed below.

• (i)

  τ_lib/2 > τ_vis (below the solid white line): f_red shows no dependence on alpha turbulent viscosity and planet mass, but it decreases with radial distance. This can be explained as follows: in this regime, viscous diffusion implies that I_v,lib ≈ I_v(r_P), and from Eq. (2) I_v,lib/I_v,cross ≈ 1 + (x_s/r_p)(3/2 – β). Eq. (17) then yields $$f_{\text{red}}\approx \frac{1}{1+4\pi \Sigma x_{\text{s}}^2\left(\frac{3}{2}-\beta \right)/M_{\text{P}}},$$(20)

  which is indeed independent of α_vis and M_P (since $x_{\text{s}}\propto \sqrt{M_{\text{P}}}$) and decreases with radial distance for our radial profiles of density and temperature.

• (ii)

  τ_lib/2 < τ_vis (above the solid white line): here, a memory effect sets in for the librating vortensity and we see that f_red now decreases with planet mass, but increases with alpha viscosity. We interpret this as a consequence of I_v,lib scaling as (1 – τ_memory/τ_mig)^3/2–β according to Eq. (6) with τ_memory ∝ τ_vis (Eq. (13)). Since τ_mig ∝ 1/M_P and ${\tau }_{\text{vis}}\propto M_{\text{P}}^2/{\alpha }_{\text{vis}}$, increasing the planet mass or decreasing the alpha viscosity leads to a higher I_v,lib/I_v,cross ratio, a stronger dynamical corotation torque, and, thus, a smaller reduction factor, f_red.

The panels clearly show that, overall, the reduction of the classical type I migration torque by the dynamical corotation torque is stronger at lower α_vis and lower β (as already seen in the previous sections), but also for higher planet masses and larger orbital distances. This is probably best seen by looking at the position of the dash-dotted curve throughout the panels, which marks where the dynamical corotation torque reduces the classical type I migration torque by a factor of 2.

Classical type I migration occurs so long as the planet wakes do not significantly perturb the disc’s background density and/or temperature profiles. The threshold mass of a planet above which its wakes start altering the background disc depends on the sound speed, the turbulent viscosity, and the relative velocity between the disc and planet (it regulates how frequently given fluid parcels are shocked through the wakes; see e.g. Wafflard-Fernandez & Baruteau 2020). Different indicators can be used to estimate this threshold mass. One is the thermal mass, which is the mass from which the planet’s Hill radius exceeds a fraction of the pressure scale height; we simply write this as M_thermal = h³M_★. Another is the feedback mass, which is given by M_feedback ≈ 3.8 × (Q/h)^–5/13 × 2/3 M_thermal (Eq. 53 in Rafikov 2002). It corresponds to the mass from which migration becomes impacted by the opening of a small gap around the planet orbit, which often comes with the formation of vortices due to the Rossby-wave instability setting at the gap edges, in particular in low-viscosity discs (McNally et al. 2019; Ziampras et al. 2024b, 2025). We note that the feedback mass expression has been derived for inviscid discs and discards the disc-planet relative drift (due e.g. to planet migration). The solid and dashed black curves in Figure 5 show where the planet mass exceeds the thermal mass and the feedback mass, respectively.

The outcome of global models of planet formation and evolution and planetary population syntheses using such models are sensitive to the planet’s growth and migration timescales (Burn & Mordasini 2024). It is a well-known problem for low-mass planets, which pose the challenge of growing faster than they migrate inward. This situation is also relevant to our study as the dynamical corotation torque comes into play once the planet has grown to several Earth masses and/or remains at fairly large orbital separations (≳ 10 au). A detailed investigation of the interplay between growth and migration via global models including our formulation of the dynamical corotation torque is the subject of a future study.

6 Discussion

6.1 On the 2D simulations

In this paper, we present the results of 2D hydrodynamical simulations for a specific physical model and numerical setup, and we now discuss how these impact our results.

Effects of grid resolution: the 2D simulations presented in Section 3 have (N_r, Ν_ϕ) = (782, 2346), with grid cells spaced linearly or logarithmically along the radial direction. The full width of the planet’s HS region is initially resolved by about 10 cells, which is sufficient to capture the effects of vortensity mixing in the HS region within a reasonable runtime. However, a finite grid resolution inevitably leads to numerical diffusion effects, which can be seen as an equivalent numerical viscosity. To see how it impacts our results, we ran additional test simulations with Q₀ = 8, β = 1 for α_vis = {0, 10^–7, 10^–6}. The simulations with α_vis = 10^–7 and α_vis = 0 show identical results that both deviate from those of the run with α_vis = 10^–6. We thus conclude that for our setup, the effects of numerical diffusion can be seen as an equivalent alpha viscosity whose magnitude is between 10^–7 and 10^–6. Furthermore, as the planet migrates inwards, its HS region gets narrower $\left(x_{\text{S}}\propto r_{\text{P}}^{1-f/2}\right)$ and less well resolved for a grid with linear radial spacing. Again, additional test simulations have shown that a grid with logarithmic radial spacing had to be used in our runs where the planet comes close to the grid’s inner wave-damping zone (which is the case of our default Q₀ = 8 runs for β = {3/2, 2}), more especially at low viscosities. This explains the grid choices made in our set of simulations (see Table 1).

Non-isothermal simulations: our 2D simulations have assessed the evolution of the disc’s vortensity in the planet’s HS region in globally isothermal discs. As already stated in Section 3.1, this choice was motivated by the need to avoid strong vortensity perturbations close to the planet’s location, and thus aimed to facilitate the (initial) estimation of I_v,lib and I_v,cross from the simulations results. However, protoplanetary discs are not thought to have a uniformly radial temperature profile, and we therefore carried out a set of additional 2D and 1D simulations with a temperature profile of T(r) ∝ r^–1 (the disc’s flaring index thus cancels out) for α_vis ∈ {10^–5, 10^–4} and β ∈ {1/2, 1, 3/2} (all 2D simulations use a grid with logarithmic radial spacing). The results of these simulations are displayed in Figure 6 and show an overall good agreement between the 1D and 2D runs, both for the I_v,lib/I_v,cross ratios and the migration rates. It indicates that our formulation of the dynamical corotation torque is not restricted to uniform temperature profiles. Still, we stress that the present study assumes discs described by a locally isothermal equation of state, in which the disc’s entropy gradient has no impact on the static and dynamical corotation torques. Future work will aim to examine how to extend our new formulation for the dynamical corotation torque to radiative disc models (discs described by an energy equation with heating and cooling source terms).

Fig. 5 Maps showing reduction factor fred, defined by Eq. (16), as a function of planet mass and orbital distance. The quantity fred quantifies how much the classical type I migration torque (sum of the Lindblad and static corotation torques) is reduced by the dynamical corotation torque: the latter is negligible where fred → 1 and strong where fred ≪ 1. Maps are obtained for four disc models that differ by the choice of the alpha turbulent viscosity (αvis) and the slope of the density profile (–β). Contour lines are plotted for fred ∈ [0.3, 0.5, 0.7, 0.9] for better visualisation. The solid white line marks the transition where τlib/2 = rvis (see text). The dashed black line marks the feedback mass (see text). Regions where the planet mass exceeds the thermal mass (equivalent to qP/h3(rP) > 1) are excluded as we expect substantial gap opening to occur in low-viscosity discs there.

6.2 An alternative 1D model for the dynamical corotation torque

As recalled in Section 2, the dynamical corotation torque depends on the inverse vortensities of the librating flow (I_v,lib) and the orbit-crossing flow (I_v,cross). Its inclusion in 1D simulations of planet formation and evolution necessitates expressions for I_v,lib and I_v,cross, which is relatively straightforward for I_v,cross (Eq. (2)), but much less straightforward for I_v,lib since the latter depends on the planet’s migration history. We show that the time evolution of I_v,lib in 2D simulations of disc-planet interactions is well reproduced by a memory-timescale argument, and the results of 2D simulations were used to estimate the memory timescale.

During the course of this study, we came up with another simple, empirical model to evaluate I_v,lib and I_v,cross, which can be expressed as follows: $${\omega }_{\text{lib}}\left(t\right)=\frac{{\tau }_{\text{vis}}{\omega }_0\left(r_{\text{P}}\left(t=0\right)\right)+t\omega \left(r_{\text{P}}\left(t\right)\right)}{{\tau }_{\text{vis}}+t},$$(21) $${\omega }_{\text{cross}}\left(t\right)=\omega \left(r_{\text{P}}\left(t\right)-x_{\text{s}}\left(t\right)\right),$$(22)

where ω₀ and ω denote the initial and instantaneous vortensity profiles, τ_vis is the viscous timescale across the planet’s HS region (see Section 2.2), and $I_{v,lib}\left(t\right)={\omega }_{\text{lib}}^{-1}\left(t\right)$ and $I_{v,\text{cross}}\left(t\right)={\omega }_{\text{cross}}^{-1}\left(t\right)$. Eq. (22) is simply the same as Eq. (2). As for Eq. (21), the rationale behind the linear combination is as follows. When t ≪ τ_vis (i.e. for very low turbulent viscosities), ω_lib(t) ≈ ω₀(r_P(t = 0)), the librating flow keeps its initial vortensity, as expected when vortensity advection dominates over vortensity diffusion. When t ≫ τ_vis, ω_lib(t) ≈ ω(r_P(t)), the vortensity of the librating flow adjusts to the instantaneous local vortensity, as expected when vortensity diffusion dominates over vortensity advection. Eqs. (21) and (22) yield the following limits for I_v,lib/I_{v, cross}: $${\frac{I_{v,\text{lib}}}{I_{v,\text{cross}}}\approx 1-\frac{x_{\text{s}}\left(t\right)}{r_{\text{P}}\left(t\right)}\frac{\text{d}\log \omega }{\text{d}\log r}|}_{r_{\text{P}}\left(t\right)}\text{for}{\tau }_{\text{vis}}\ll t,$$(23) $$\approx \frac{\omega \left(r_{\text{P}}\left(t\right)-x_{\text{s}}\left(t\right)\right)}{{\omega }_0\left(r_{\text{P}}\left(t=0\right)\right)}\text{for}{\tau }_{\text{vis}}\gg t.$$(24)

In the first limit (Eq. (23)), I_v,lib/I_v,cross is entirely set by the background vortensity gradient, while the second limit (Eq. (24)) corresponds to the case where the librating flow has infinite vortensity memory.

Figure 7 compares the results of the 2D simulations presented in Section 3.2 with those of 1D simulations using the same setup as in Section 4.1, but with I_v,lib and I_v,cross computed from Eqs. (21) and (22). We see that the alternative model proposed in this section also works very well overall. We see that for β < 3/2, the alternative model better reproduces the initial increase in I_v,lib/I_v,cross but slightly overestimates its stationary value, thereby predicting slightly slower inward migration than with our fiducial memory model. The 1D simulations for β > 3/2 show quite remarkable agreement with the 2D simulations.

6.3 Comparison to previous work

A simple model for evaluating I_v,lib/I_v,cross was derived in Paardekooper (2014) based on the evolution towards a steady state of the vortensity in the planet’s HS region due to viscous diffusion and planet migration (see Eq. (28) of Paardekooper 2014): $$1-\frac{I_{v,\text{lib}}}{I_{v,\text{cross}}}=\left(\frac{3}{2}-\beta \right)\frac{x_{\text{s}}^2}{6r_{\text{P}}v}{\stackrel{·}{r}}_{\text{P}},$$(25)

which can be recast as $$\frac{I_{v,\text{lib}}}{I_{v,\text{cross}}}=1-\left(\frac{3}{2}-\beta \right)\frac{{\tau }_{\text{vis}}}{6{\tau }_{\text{mig}}}.$$(26)

When τ_lib ≪ τ_vis ≪ |τ_mig|, the memory timescale in Eq. (13) is ≈0.4τ_vis = 2τ_vis/5, and Eq. (6) yields $$\frac{I_{v,\text{lib}}}{I_{v,\text{cross}}}\approx 1-\left(\frac{3}{2}-\beta \right)\frac{2{\tau }_{\text{vis}}}{5{\tau }_{\text{mig}}},$$(27)

which is similar to the original model of Paardekooper (2014), except for the constant factor in front of τ_vis/τ_mig. We tested this model by rerunning 1D simulations with the same setup as in Section 4.1, using the expression in Eq. (14) for the dynamical corotation torque, which is valid when τ_lib ≪ |τ_mig|, and using the expression in Eq. (26) for I_v,lib/I_v,cross.

Figure 8 compares the results of these 1D simulations with those of the 2D simulations presented in Section 3.2. We see that the 1D simulations underestimate the ratio I_v,lib/I_v,cross and therefore predict faster inward migration than observed in our 2D simulations. This may be related, at least partly, to the constant factor in front of τ_vis/τ_mig (1/6 vs. 2/5 between Eqs. (26) and (27)). Also, since Eq. (26) was derived as a stationary solution, it cannot reproduce the gradual increase in I_v,lib/I_v,cross before reaching a quasi-steady state. Finally, we see that the results of the 1D runs significantly depart from those of the 2D runs for β > 3/2. This is not surprising since the above formulation for I_v,lib/I_v,cross assumes a steady migration rate, which does not apply to the runaway migration regime. We briefly note that, due to a fast increase in the migration rate for α_vis = 10^–5 and β = 2, our 1D simulation showed unstable behaviour and was stopped after ~500 orbits. To our knowledge, the above results constitute the first systematic benchmark of the dynamical corotation torque model of Paardekooper (2014) against hydrodynamical simulations of disc-planet interactions.

Since the study of Paardekooper (2014), the growth of interest in understanding planet migration and evolution in low-viscosity discs has brought up a number of works on the dynamical corotation torque via 2D and 3D (magneto) hydro-dynamical simulations and simple 1D simulations. Here, we briefly summarise their results. A first attempt at including the effects of the dynamical corotation torque in 1D planet formation models was done by Ogihara et al. (2017) to study the effects of magnetised disc winds on type I migration. They used the dynamical corotation torque expression of Paardekooper (2014) with I_v,lib/I_v,cross given by Eq. (26), but where τ_mig is the characteristic timescale for the relative migration between the planet and disc. They predicted that super-Earths would have a bimodal distribution of semi-major axes due to how winds affect the disc structure and the resulting type I migration. Ndugu et al. (2021) investigated the dynamical corotation torque by assuming conservation of vortensity in the horseshoe region but calculating the evolution of the dynamical corotation torque by evolving the mass of the horseshoe region (including effects of gap opening, viscous spreading and planet accretion). These authors concluded that planet accretion from the horseshoe region and the dynamical corotation torque can efficiently slow down type I migration. In McNally et al. (2017) and McNally et al. (2018), the effect of a net radial in- or outflow induced by Maxwell stresses in low-viscosity discs was also investigated. They found a dynamical corotation torque arising from the flow-induced asymmetric distortion of the corotation region, leading to four distinct migration regimes depending on the setup (i.e. locked migration to the inward or outward radial gas flow and outward runaway migration). They found that runaway outward migration saturates to a steady speed and provided an analytical model. In McNally et al. (2020), it was shown that a wind-driven laminar accretion flow through the disc surface layer does not significantly impact the torques exerted onto embedded planets. However, they also showed that the dynamical corotation torque can change behaviour dramatically when going from 2D to 3D, switching signs and leading to fast inward migration due to the dissipation of buoyancy waves (see also Zhu et al. 2012).

Further 3D inviscid radiation hydrodynamic simulations showed that radiative cooling and diffusion can damp buoyancy waves more or less efficiently depending on the disc’s local physical properties (Yun et al. 2022; Ziampras et al. 2024a), which affects the sign and magnitude of the dynamical corotation torque. All the aforementioned works show that a complete picture of the dynamical corotation torque has not yet emerged and that there is still much room for the improvement of its formulation for global models of planet formation and evolution, which can serve, via population syntheses, to compare the different models with observations.

Fig. 6 Same as Figure 2, but for additional simulations with a background temperature profile decreasing as 1/r, demonstrating the ability of our simple 1D model to also reproduce non-isothermal simulation results.

Fig. 7 Same as Figure 2, except that 1D simulations use alternative formulation for dynamical corotation torque of Section 6.2. It shows remarkable agreement for the runaway migration (β > 3/2) and good overall agreement for the slow-down migration simulations. For an easier direct comparison of the most relevant simulations, see right column in Figure B.1.

Fig. 8 Same as Figure 2, except that 1D simulations use original formulation for dynamical corotation torque of Paardekooper (2014), which is recalled in Section 6.3. Note that this steady-state approach is not able to reproduce the initial build-up of vortensity ratio and the final values are generally under-predicted.

7 Summary and conclusions

The dynamical corotation torque can greatly impact the migration of low-mass planets in their protoplanetary disc, particularly for low-viscosity wind-driven discs. Yet, it is often discarded in global models of planet formation and evolution, which rely on 1D simulations. Its calculation indeed requires analytical prescriptions for the vortensity of the librating and orbit-crossing flows in the planet’s horseshoe region, which, for the librating flow, depend on the migration history of the planet. In this work, we proposed two simple prescriptions for the time evolution of the librating vortensity.

Our main prescription is based on a memory effect (Section 2.2), which features a memory timescale on which the librating flow forgets about its initial vortensity and adjusts to that of the ambient gas, after which the ratio of inverse vorten-sities between the librating and orbit-crossing flows becomes stationary. 2D hydrodynamical simulations assuming locally isothermal discs were used to estimate this memory timescale, which is a fraction of the viscous timescale across the horseshoe region. An alternative prescription for the inverse vortensities of the librating and orbit-crossing flows is proposed in Section 6.2.

One-dimensional simulations including the aforementioned prescriptions were carried out and compared to 2D hydrody-namical simulations for locally isothermal discs. The planet’s orbital evolutions in the 1D and 2D simulations are in very good agreement for a wide parameter space, whether the dynamical corotation torque slows down inward migration or promotes runaway inward migration. For typical shallow density and temperature profiles in discs, we provide maps showing how much the dynamical corotation torque reduces the classical type I migration torque (sum of the Lindblad and static corotation torques) as a function of planet mass and orbital distance (Fig. 5, Section 5). This reduction is larger at higher planet mass and larger orbital distance. For a young disc with a surface density profile in r^–1/2 and α_vis = 10^–4, it is by a factor ≈50% for a 10-Earth-mass planet at 10 au.

In summary, the dynamical corotation torque should be taken into account in global models of planet formation and evolution for discs with low turbulent viscosity (α_vis ≲ 10^–4). The presented framework of using hydrodynamical simulations of disc-planet interactions to inform 1D analytical prescriptions for the vortensity in the planet’s horseshoe region will in the future help assess the impact of other physical processes on the dynamical corotation torque such as vortensity production due to baroclinic forcing (Ziampras et al. 2024a,b).

Acknowledgements

We thank Oliver Schib, Richard Nelson, Aurélien Crida, Alessandro Morbidelli, Nicolas Kaufmann, and Andrin Kessler for stimulating discussions on this work. We further thank Alessandro Morbidelli for a thorough reading of a first version of this manuscript. We also thank the anonymous referee for their comments that helped improve the manuscript. J.W. and C.M. acknowledge the support from the Swiss National Science Foundation under grant 200021_204847 “PlanetsInTime”. Part of this work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. Calculations were performed on the Horus cluster of the Division of Space Research and Planetary Sciences at the University of Bern. Part of the numerical simulations were performed on the CALMIP Supercomputing Centre of the University of Toulouse.

Appendix A Additional simulations with higher aspect ratio

In the main text, we present results of simulations for inwardly migrating low-mass planets starting at a disc’s aspect ratio h₀ = 0.05. Since our fiducial disc model assumes h(r) ∝ r^1/2, our simulations automatically tested the impact of having a smaller aspect ratio through migration (down to h(0.5r₀) ~ 3.5%). In order to further test our 1D models of the dynamical corotation torque against larger aspect ratios, we have carried out additional simulations for h₀ = 0.07. We have increased the planet mass to maintain the same $q_{\text{p}}/h_0^3$ ratio as in our fiducial model, such that the density perturbations associated to the planet wakes reach very similar values. We have also increased Σ₀ to retain Q₀ = 8. The initial migration rate, which scales as q_p/{h₀Q₀], is about twice as large as in our fiducial model for the same Q₀. These additional simulations used α_vis = 10^–5 and $\beta =\left\{\frac{1}{2},1\right\}$. Results are shown in Figure A.1.

The memory timescale model (Section 4.1) predicts a slightly different time evolution for the inverse vortensity ratio (see bottom-left panel), but the magnitude is well reproduced and therefore so is the migration rate. The alternative model of Section 6.2 reproduces nearly perfectly the time evolution of the inverse vortensity ratio, the migration rate and the planet’s orbital distance. The steady-state solution model discussed in Section 6.3 underpredicts the inverse vortensity ratio, leading to too high migration rates.

We finally note that these simulations all show a constantly increasing inverse vortensity ratio. In both our memory timescale model and the steady-state solution model, I_v,lib/I_v,cross – 1 scales with τ_vis/τ_mig (see Eqs. 26 and 27; these are valid for τ_lib ≪ τ_vis ≪ |τ_mig|, and the latter inequality is marginally fulfilled in the simulations presented here). Since ${\tau }_{\text{vis}}\equiv q_{\text{p}}/\left\{\alpha h^3\left(r_{\text{p}}\right)\right\}\times {\Omega }^{-1}\propto r_{\text{p}}^{3/2-3f},{\tau }_{\text{vis}}$ remains stationary for f = 1/2. When considering only the Lindblad torque for simplicity, |τ_mig| ∝ r^1/2+β for f = 1/2: |τ_mig| decreases as the planet migrates inwards. This explains why I_v,lib/I_v,cross also increases with time. Similarly, in the alternative model of Section 6.2, so long as t ≪ τ_vis, I_v,lib/I_v,cross is given by Eq. (24) and smoothly increases with time as the planet migrates inwards. It will change behaviour only if t ≫ τ_vis. At the end of the simulations, τ_vis/t does not go below ≈ 0.5. The final value for I_v,lib/I_v,cross is larger in the h₀ = 0.07 simulations than in the h₀ = 0.05 runs because the planet has migrated farther inwards.

Fig. A.1 Direct comparison of 2D and 1D simulations results for h0 = 0.07, qp = 2.7 × 10–5 and Q0 = 8. Rows are the same as in Fig. 2. Columns correspond to different approaches for calculating the inverse vortensities Iv,lib and Iv,cross: the memory timescale model (Section 4.1, left), the alternative model (Section 6.2, middle) and the steady state solution (Section 6.3, right). Opaque lines are for 1D simulations, transparent lines for 2D simulations.

Appendix B Alternative figure comparing results of 2D and 1D simulations

Throughout the manuscript, results of 2D and 1D simulations are compared for different 1D models of the dynamical corotation torque. A side-by-side comparison between 2D and 1D results is displayed in Figs. 2, 6, 7 and 8. These figures also help assess the impact of including or not the dynamical corotation torque in 1D models. This appendix is used to include Fig. B.1, which provides a more direct comparison between the results of 2D and 1D simulations by overlying them on same panels. The left column shows the results of the memory timescale model in Fig. 2, the right column those of the alternative model in Fig. 7.

Fig. B.1 Overlay of 2D and 1D simulations results from Fig. 2 (memory timescale model of Section 4.1, left column) and Fig. 7 (alternative model of Section 6.2, right column). Rows are the same as in Fig. 2. Opaque lines display 1D simulations results, whereas transparent lines show the results of 2D hydrodynamical simulations. Cases with αvis = 10–3 are only shown for β ≥ 1.5 to make plots less crowded.

References

All Tables

All Figures

	Fig. 1 Sketch illustrating perturbed vortensity around planet’s vicinity (black circle) with relevant flows and torques highlighted: the inner and outer spiral density waves inducing the Lindblad torque, horseshoe U-turns ahead of and behind the planet leading to the static corotation torque, and the librating and orbit-crossing flows giving rise to the dynamical corotation torque.
In the text

Fig. 2 From top to bottom: time evolution of planet’s migration rate, orbital distance, and inverse vortensity ratio Iv,lib/Iv,cross for simulations with Q0 = 8. The results of our 2D (with FARGO3D) and 1D (with the Bern model) simulations are shown side by side for easy direct comparison. Lines are coloured by αvis and different line styles indicate varying slopes of the surface density, β. The migration rates of the 2D simulations (top left panel) are shown as moving averages over 100 orbits. For the 1D simulations (right), we additionally show results of simulations that discard the dynamical corotation torque as transparent lines. For an easier direct comparison of the most relevant simulations, see the left column in Figure B.1 of Appendix B.

In the text

Fig. 3 Evolution of disc vortensity for simulation with Q0 = 8, β = 1 and αvis = 10–5. The top panels show snapshots of the vortensity field around the planet at 500, 1500, and 2500 orbits. The vortensity deficit in the HS region, and the vortensity advection of the orbit-crossing flow from the inner to the outer edges of the HS region are clearly visible. The vortensities of the orbit-crossing and librating flows measured by the method described in Section 3.2.2 are shown as contour lines in pink and cyan, respectively. The grey vertical lines mark the orbital radius of the planet and the extent of its horseshoe region (see Section 2.1 for the expression of the horseshoe half-width xs). The bottom panel shows the time evolution of the measured vortensities for the orbit-crossing and librating flows. Crosses mark the results shown in the top panels. The time evolution of Iv,lib/Iv,cross is shown in black.

In the text

	Fig. 4 Left and middle panels: time evolution of Iv,lib/Iv,cross in simulations with varying Q0 ∈ {8, 12, 16, 24, 32}, $\beta \in \left\{\frac{1}{2},1\right\}$, for αvis = 10–4 (left) and αvis = 10–5 (middle). Right panel: time evolution of ratio between memory timescale and viscous timescale across HS region for αvis = 10–5 simulations.
In the text

Fig. 5 Maps showing reduction factor fred, defined by Eq. (16), as a function of planet mass and orbital distance. The quantity fred quantifies how much the classical type I migration torque (sum of the Lindblad and static corotation torques) is reduced by the dynamical corotation torque: the latter is negligible where fred → 1 and strong where fred ≪ 1. Maps are obtained for four disc models that differ by the choice of the alpha turbulent viscosity (αvis) and the slope of the density profile (–β). Contour lines are plotted for fred ∈ [0.3, 0.5, 0.7, 0.9] for better visualisation. The solid white line marks the transition where τlib/2 = rvis (see text). The dashed black line marks the feedback mass (see text). Regions where the planet mass exceeds the thermal mass (equivalent to qP/h3(rP) > 1) are excluded as we expect substantial gap opening to occur in low-viscosity discs there.

In the text

	Fig. 6 Same as Figure 2, but for additional simulations with a background temperature profile decreasing as 1/r, demonstrating the ability of our simple 1D model to also reproduce non-isothermal simulation results.
In the text

	Fig. 7 Same as Figure 2, except that 1D simulations use alternative formulation for dynamical corotation torque of Section 6.2. It shows remarkable agreement for the runaway migration (β > 3/2) and good overall agreement for the slow-down migration simulations. For an easier direct comparison of the most relevant simulations, see right column in Figure B.1.
In the text

	Fig. 8 Same as Figure 2, except that 1D simulations use original formulation for dynamical corotation torque of Paardekooper (2014), which is recalled in Section 6.3. Note that this steady-state approach is not able to reproduce the initial build-up of vortensity ratio and the final values are generally under-predicted.
In the text

Fig. A.1 Direct comparison of 2D and 1D simulations results for h0 = 0.07, qp = 2.7 × 10–5 and Q0 = 8. Rows are the same as in Fig. 2. Columns correspond to different approaches for calculating the inverse vortensities Iv,lib and Iv,cross: the memory timescale model (Section 4.1, left), the alternative model (Section 6.2, middle) and the steady state solution (Section 6.3, right). Opaque lines are for 1D simulations, transparent lines for 2D simulations.

In the text

Fig. B.1 Overlay of 2D and 1D simulations results from Fig. 2 (memory timescale model of Section 4.1, left column) and Fig. 7 (alternative model of Section 6.2, right column). Rows are the same as in Fig. 2. Opaque lines display 1D simulations results, whereas transparent lines show the results of 2D hydrodynamical simulations. Cases with αvis = 10–3 are only shown for β ≥ 1.5 to make plots less crowded.

In the text