# The importance of the dynamical corotation torque for the migration of low-mass planets

Jesse Weder, Clément Baruteau, Christoph Mordasini

► To cite this version:

Jesse Weder, Clément Baruteau, Christoph Mordasini. The importance of the dynamical corotation torque for the migration of low-mass planets. Astronomy & Astrophysics - A&A, 2025, 701, pp.A16. 10.1051/0004-6361/202554620 . hal-05401304

HAL Id: hal-05401304

https://hal.science/hal-05401304v1

Submitted on 20 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# The importance of the dynamical corotation torque for the migration of low-mass planets

## 1D analytical prescriptions verified by 2D hydrodynamical simulations

Jesse Weder<fnref n="1" />,*, Clément Baruteau<fnref n="2" />, and Christoph Mordasini<fnref n="1" />

Received 18 March 2025 / Accepted 27 June 2025

## ABSTRACT

*Context.* Recent developments in the evolution of protoplanetary discs have suggested that planet formation occurs in regions of the discs with low levels of turbulent viscosity. In such environments, the dynamical corotation torque is thought to play an important role by slowing down the migration of low-mass planets (type I migration).

*Aims.* We aim to provide a simple analytical prescription for the dynamical corotation torque for use in 1D global models of planet formation and evolution and assess the importance of the dynamical corotation torque for the migration of low-mass planets in low-viscosity discs.

*Methods.* We propose simple prescriptions for calculating in 1D the time evolution of the vortilities of the librating and orbit-crossing flows around a low-mass planet, which both enter the analytical expression for the dynamical corotation torque. One of our prescriptions involves a memory timescale for the librating flow, and 2D hydrodynamical simulations of disc-planet interactions were used to assess the memory timescale and validate our model.

*Results.* The orbital evolution of a low-mass planet is calculated using 1D simulations, where the dynamical corotation torque features our prescriptions for the vortilities of the librating and orbit-crossing flows, and using 2D hydrodynamical simulations of disc-planet interactions, assuming locally isothermal discs. We find very good agreement between the 1D and 2D simulations for a wide parameter space, whether the dynamical corotation torque slows down or accelerates inward migration. We provide maps showing how much the dynamical corotation torque reduces the classical type I migration torque as a function of planet mass and orbital distance. The reduction is about 50% for a 10-Earth-mass planet at 10 au in a young disc with a surface density profile in $r^{-1/2}$ and an alpha viscosity of $10^{-4}$.

*Conclusions.* In discs with low turbulent viscosity, the dynamical corotation torque should be taken into account in global models of planet formation and evolution as it can strongly slow down type I migration.

**Key words.** hydrodynamics – planets and satellites: formation – protoplanetary disks – planet-disk interactions

## 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 $\alpha_{vis} \gtrsim$ a few $\times 10^{-3}$ (Bitsch et al. 2015; Liu et al. 2019; Emsenhuber et al. 2021; Dąż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 $\alpha_{vis} \lesssim 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 $\alpha_{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 magnetohydrodynamical (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 $\alpha_{vis} \lesssim 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 $\alpha_{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

<footnote n="1">¹ Division of Space Research and Planetary Sciences, Physics Institute, University of Bern, Gesellschaftsstrasse 6, 3012 Bern, Switzerland</footnote>

<footnote n="2">² IRAP, Université de Toulouse, CNRS, Université Paul Sabatier, CNES, Toulouse, France</footnote>

<footnote>* Corresponding author; jesse.weder@unibe.ch</footnote>
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 parameterisations 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<fnref n="1" /> the Lindblad torque ($\Gamma_{\text{lin}}$, angular momentum exchange between the planet and its wakes in the disc) and the static corotation torque ($\Gamma_{\text{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 ($\Gamma_{\text{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.

<figure><img src="image_2.png" /><figcaption>**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.</figcaption></figure>

Following Paardekooper (2014) and McNally et al. (2017), the dynamical corotation torque can be written as

$$ \Gamma_{\text{c,dyn}} = \frac{3}{4} r_P \Omega_P^3 \int_{y_s}^{x_s} (\mathcal{I}_{\text{v,lib}} - \mathcal{I}_{\text{v,cross}}) x^2 dx, \quad (1) $$

where $r_P$ and $\Omega_P$ denote the planet's orbital radius and itsKeplerian frequency, and $\mathcal{I}_{\text{v,lib}}$ and $\mathcal{I}_{\text{v,cross}}$ give the inverse vortensity of the librating and orbit-crossing flows. The quantity $\mathcal{I}_v$ is defined as $\mathcal{I}_v = \Sigma/\{\nabla \times \mathbf{v}\} \cdot \mathbf{e}_z\}$, with $\Sigma$ the surface density of the disc, $\mathbf{v}$ its velocity field, and $\mathbf{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 + \dot{r}_P\tau_{\text{lib}}/2$, with $\dot{r}_P$ being the planet's migration rate and $\tau_{\text{lib}} = 8\pi r_P/(3\Omega_P 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.

### 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 \ll h^3(r_P)$ –where h is the disc's

<footnote n="1">¹ Thermal torques that arise from accretion onto the planet (Masset 2017) and dust-related torques (Benítez-Llambay & Pessah 2018; Guilera et al. 2023; Chrenko et al. 2024) are discarded here for simplicity.</footnote>
aspect ratio (Lin & Papaloizou 1993)– the inverse vortensity of the orbit-crossing flow takes its local, unperturbed value. This results in

$$ \mathcal{I}_{v,\text{cross}} = \mathcal{I}_v(r_p - x_s) \approx \mathcal{I}_v(r_p) \times \left[ 1 - \frac{x_s}{r_p} \left( \frac{3}{2} - \beta \right) \right], \quad (2) $$

where we assumed that the disc’s surface density $\Sigma \propto r^{-\beta}$ and that the disc’s velocity is purely Keplerian<fnref n="2" />. The half-width of the planet’s horseshoe region is computed as $x_s = Cr_p \sqrt{q_p/h(r_p)}$, 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 $(\nabla\Sigma \times \nabla P)/\Sigma^3$ (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 \times \partial\Sigma/\partial\varphi$ 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, $\tau_{\text{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 - \tau_{\text{memory}}$, when it was located at a radial distance of $r_p - \dot{r}_p\tau_{\text{memory}}$. More specifically, we assume that

$$ \mathcal{I}_{v,\text{lib}} = \mathcal{I}_v(\min[r_p - \dot{r}_p\tau_{\text{memory}}, r_0]), \quad (3) $$

where the minimum function $\min$ ensures that the vortensity cannot originate from further away than the planet’s initial location $r_0$. As mentioned already, in this work we limited ourselves to inward-migrating planets ($\dot{r}_p < 0$), and for outward migration ($\dot{r}_p > 0$) one should replace the minimum function in Eq. (3) with the maximum one. Introducing the migration timescale $\tau_{\text{mig}}$ such that $\dot{r}_p = r_p/\tau_{\text{mig}}$, one obtains

$$ \mathcal{I}_{v,\text{lib}} = \mathcal{I}_v(\min[r_p \times (1 - \tau_{\text{memory}}/\tau_{\text{mig}}), r_0]) \quad (4) $$

$$ = \mathcal{I}_v \left( r_p \times \min \left[ 1 - \frac{\tau_{\text{memory}}}{\tau_{\text{mig}}}, \frac{r_0}{r_p} \right] \right) \quad (5) $$

$$ = \mathcal{I}_v(r_p) \times \left( \min \left[ 1 - \frac{\tau_{\text{memory}}}{\tau_{\text{mig}}}, \frac{r_0}{r_p} \right] \right)^{3/2-\beta} . \quad (6) $$

Viscous diffusion will take $\mathcal{I}_{v,\text{lib}}$ towards the local unperturbed vortensity $\mathcal{I}_v(r_p)$ on a timescale similar to the viscous timescale over the HS region $\tau_{\text{vis}} = x_s^2/v$. For high turbulent viscosity $v$, the viscous timescale is short compared to the migration timescale, and $\mathcal{I}_{v,\text{lib}}$ is close to the local unperturbed value $\mathcal{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 $\mathcal{I}_{v,\text{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 ($\beta = 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 $\tau_{\text{memory}}$. We expect $\tau_{\text{vis}}$ to pose an upper limit on $\tau_{\text{memory}}$. However, this needs to be confirmed by measuring the vortensity evolution in 2D hydrodynamical 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, \phi]$ with $r \in [0.4, 2.5]r_0$ and $\phi \in [-\pi, \pi]$, where $r_0$, 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 $N_\phi = 2346$ cells linearly spaced along $\phi$. 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 \in [0.4, 0.52] \cup [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, $\Sigma_0(r) = \Sigma_0(r/r_0)^{-\beta}$, with $\Sigma_0$ the initial surface density at the planet’s initial location $r_0$. We assume the disc evolution is driven by radial transport of angular momentum and is parameterised by the classical $\alpha$-viscosity approach (Shakura & Sunyaev (1973); the alpha viscosity is taken uniform and denoted by $\alpha_{\text{vis}}$). The disc’s aspect ratio is given by $h(r) = h_0(r/r_0)^f$ with $h_0 = h(r_0)$. In our fiducial simulations we adopted $h_0 = 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 $\propto \partial_\phi \Sigma \partial_r T$ 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 = 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_\star \approx 10^{-5}$, which corresponds to $M_p \approx 3.3 M_\oplus$ ($M_\oplus$ 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^3 \ll 1$. We used a

<footnote n="2">² The vortensity of a Keplerian disc is $(\nabla \times \mathbf{v}) \cdot \mathbf{e}_s / \Sigma = \frac{1}{r_\Sigma} \frac{\partial}{\partial r} (r^2 \Omega) = \frac{\Omega}{2\Sigma}$.</footnote>
<table><caption>**Table 1.** Summary of initial parameters explored in 2D hydrodynamical simulations presented in Section 3.2 (Q₀: Toomre-Q parameter at the planet's initial orbital radius; α<sub>vis</sub>: turbulent alpha viscosity, β = -d log Σ₀/d log r).</caption><thead><tr><th>Q<sub>0</sub></th><th>α<sub>vis</sub></th><th>β</th><th>Grid's radial spacing</th></tr></thead><tbody><tr><td>8</td><td>{10<sup>-3</sup>, 10<sup>-4</sup>, 10<sup>-5</sup>}</td><td>{1/2, 1}</td><td>Linear</td></tr><tr><td></td><td>{10<sup>-3</sup>, 10<sup>-4</sup>, 10<sup>-5</sup>}</td><td>{3/2, 2}</td><td>Logarithmic</td></tr><tr><td>12</td><td>{10<sup>-4</sup>, 10<sup>-5</sup>}</td><td>{1/2, 1}</td><td>Linear</td></tr><tr><td>16</td><td>{10<sup>-4</sup>, 10<sup>-5</sup>}</td><td>{1/2, 1}</td><td>Linear</td></tr><tr><td>24</td><td>{10<sup>-4</sup>, 10<sup>-5</sup>}</td><td>{1/2, 1}</td><td>Linear</td></tr><tr><td>32</td><td>{10<sup>-4</sup>, 10<sup>-5</sup>}</td><td>{1/2, 1}</td><td>Linear</td></tr></tbody></table>

smoothing length of 0.6H(r<sub>p</sub>) 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 α<sub>vis</sub>, β, 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(r_0)}{\pi G \Sigma_0} \qquad (7)$$

Our default value for Q₀ is 8, for which we ran simulations with α<sub>vis</sub> ∈ {10<sup>-3</sup>, 10<sup>-4</sup>, 10<sup>-5</sup>} and β ∈ {1/2, 1, 3/2, 2}. 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}, α<sub>vis</sub> ∈ {10<sup>-4</sup>, 10<sup>-5</sup>}, and β ∈ {1/2, 1} 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<sub>p</sub> ≲ 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 Ω⁻¹(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 α<sub>vis</sub> 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 α<sub>vis</sub> from 10⁻³ to 10⁻⁴ leads to slightly higher migration rates, and then going further down to α<sub>vis</sub> = 10⁻⁵ 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 Lindblad torque is independent of viscosity (Muto & Inutsuka 2009). Only at low α<sub>vis</sub> 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<sub>v,lib</sub> - I<sub>v,cross</sub> as it is what sets the dynamical corotation torque (see Eq. (1)). To measure I<sub>v,lib</sub> and I<sub>v,cross</sub> 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<sub>v,cross</sub> 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 φ - φₚ ∈ [-0.2, -0.1] at r = rₚ (φₚ 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<sub>v,lib</sub> is averaged over azimuth at the planet's orbital radius, excluding a small region around the planet; in practice, it was averaged over φ - φₚ ∈ [-π, -0.4] ∪ [0.4, π] at r = rₚ. 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<sub>v,lib</sub> was averaged over φ - φₚ ∈ [0.4, 0.6] at r = rₚ.

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ₚ - xₐ, which is in agreement with Eq. (2). In the lower panel of the figure,
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** From top to bottom: time evolution of planet's migration rate, orbital distance, and inverse vortensity ratio *I*<sub>v,lib</sub>/*I*<sub>v,cross</sub> for simulations with *Q*<sub>0</sub> = 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 α<sub>vis</sub>, 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.</figcaption></figure>

we stress that, although both vortivities increase with time, the
ratio *I*<sub>v,lib</sub>/*I*<sub>v,cross</sub> 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 evolu-
tion of *I*<sub>v,lib</sub>/*I*<sub>v,cross</sub> in the simulations. For *β* < 3/2, the planet

moves up the disc’s vortensity gradient and thus has a vorten-
sity deficit in its HS region, which explains why *I*<sub>v,lib</sub>/*I*<sub>v,cross</sub> >
1. Quite remarkably, in all our simulations with *β* < 3/2,
*I*<sub>v,lib</sub>/*I*<sub>v,cross</sub> reaches a near-stationary value, while the planets
keeps migrating, after a characteristic timescale that depends
on *β* and *α*<sub>vis</sub>. The lower the viscosity, the longer it takes to
reach a steady state, the higher the final value for *I*<sub>v,lib</sub>/*I*<sub>v,cross</sub>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Evolution of disc vortensity for simulation with $Q_0 = 8, \beta = 1$ and $\alpha_{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 $x_s$). 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 $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ is shown in black.</figcaption></figure>

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 $\alpha_{vis} \le 10^{-4}$. For $\beta > 3/2$, the planet moves down the vortensity gradient and thus builds up a vortensity excess in its HS region, implying that $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross} < 1$. This time, the runaway implies that $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ keeps decreasing with time and does not reach a steady state value. Finally, for $\beta = 3/2$ we see that $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ stays at unity regardless of turbulent viscosity, as expected.

### 3.3. Estimating the memory timescale

As shown in Section 3.2.2, $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ eventually reaches a quasi-steady state when migration does not run away ($\beta < 3/2$). While for simulations with $\alpha_{vis} \ge 10^{-4}$ the final value of $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ can essentially be attributed to the background

vortensity gradient, for simulations at lower viscosity ($\alpha_{vis} = 10^{-5}$) $\mathcal{I}_{v,lib}/\mathcal{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 - \tau_{memory}/\tau_{mig} < r_0/r_P$ for $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ to reach a constant value. Reversing Eq. (6) then gives

$$ \tau_{memory} = \tau_{mig} \times \left( 1 - \left( \frac{\mathcal{I}_{v,lib}}{\mathcal{I}_{v}(r_P)} \right)^{\frac{1}{3/2-\beta}} \right), \quad (8) $$

where we recall that the migration timescale $\tau_{mig} = r_P/\dot{r}_P$. Although migration rates are close to stationary in simulations with $\beta < 3/2$ (see Fig. 2), the migration timescale varies over time due to the planet moving inwards, and so does the memory timescale.
<figure><img src="image_5.png" /><figcaption>Fig. 4. Left and middle panels: time evolution of $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ in simulations with varying $Q_0 \in \{8, 12, 16, 24, 32\}$, $\beta \in \{\frac{1}{2}, 1\}$, for $\alpha_{vis} = 10^{-4}$ (left) and $\alpha_{vis} = 10^{-5}$ (middle). Right panel: time evolution of ratio between memory timescale and viscous timescale across HS region for $\alpha_{vis} = 10^{-5}$ simulations.</figcaption></figure>

We computed the memory timescale in an extended set of 2D simulations that include higher values for $Q_0$ (12, 16, 24 and 32), for $\beta \in \{\frac{1}{2}, 1\}$, and for $\alpha_{vis} \in \{10^{-4}, 10^{-5}\}$. Results are shown in Figure 4. The left panel first shows the time evolution of $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ for $\alpha_{vis} = 10^{-4}$. Two trends emerge: (i) for a given $Q_0$, $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ is larger for smaller $\beta$; and (ii) for a given $\beta$, $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ decreases with increasing $Q_0$. The viscous and libration timescales over the HS region happen to have similar values within about 30% ($\tau_{vis} \approx 130$ orbital periods at $r_0$, $\tau_{lib} \approx 100$ orbital periods at $r_0$). 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 $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ is essentially explained by the background vortensity gradient (see Eq. (2) with $\mathcal{I}_{v,lib} \sim \mathcal{I}_v(r_p)$). Evaluating Eq. (2) at the planet's initial location actually gives $\mathcal{I}_v(r_p)/\mathcal{I}_{v,cross} \approx 1.014$ for $\beta = 1/2$ and $\approx 1.007$ for $\beta = 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 $\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ for $\alpha_{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 \le 24$). Having previously shown that the background vortensity gradient also contributes to the final value for $\mathcal{I}_{v,lib}/\mathcal{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 $\beta$ values reach a quasi-steady state for different values of $\mathcal{I}_{v,lib}/\mathcal{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_0$), but overall $\tau_{memory} \sim [0.3-0.4] \times \tau_{vis}$ is a good match to our simulations results. These experiments allow us to conclude that two effects impact the ratio $\mathcal{I}_{v,lib}/\mathcal{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 $\approx 0.4\tau_{vis}$.

## 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} (r^2 \Sigma \alpha_{vis} c_s^2) \right], \quad (9) $$

with $c_s = hr\Omega$ being the sound speed and $\alpha_{vis}$ the alpha turbulent viscosity. As in our 2D simulations, the disc's aspect ratio was chosen $h = h_0(r/r_0)^f$ with $f = 1/2$ (uniform temperature profile), the initial surface density was $\Sigma_0(r/r_0)^{-\beta}$, and $\alpha_{vis}$ was taken constant. The 1D grid used 1244 cells logarithmically spaced over $r \in [0.2, 2.5]r_0$. 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^{-4}r_P$ per time step.

At the beginning of the simulations, a planet was inserted at $r_0$ 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)<fnref n="3" />. For the dynamical corotation torque, Eq. (1) gives

$$ \Gamma_{c,dyn} = \frac{1}{4} r_P \Omega_P^3 (\mathcal{I}_{v,lib} - \mathcal{I}_{v,cross}) (x_s^3 - y_s^3), \quad (10) $$

assuming $\mathcal{I}_{v,lib}$ and $\mathcal{I}_{v,cross}$ are uniform over the HS region. We further set $y_s$ such that

$$ y_s = \max(x_s + \dot{r}_P \tau_{lib}/2, 0), \quad (11) $$

which ensures that in the fast migrating regime (runaway), the dynamical torque reaches

$$ \Gamma_{c,dyn} = \frac{1}{4} r_P \Omega_P^3 (\mathcal{I}_{v,lib} - \mathcal{I}_{v,cross}) x_s^3 \quad (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 $\beta$. Finally, $\mathcal{I}_{v,cross}$ was calculated using Eq. (2) and $\mathcal{I}_{v,lib}$ via Eq. (6) with

$$ \tau_{memory} = 0.4 \tau_{vis} \times \max \left( 1 - \frac{\tau_{lib}}{2\tau_{vis}}, 0 \right). \quad (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 $\tau_{vis} \approx 2\tau_{lib}$, corresponding to $\alpha_{vis} \approx 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 $\tau_{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_0 = 8$ and the same range of values for $\beta$ (1/2,1) and $\alpha_{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 ($\beta = 3/2$), the disc torque is entirely given by the Lindblad torque and does not depend on $\alpha_{vis}$. We note that the 2D torque starts to slightly exceed the 1D torque over long timescales ($\ge 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 ($\beta < 3/2$), it is interesting to notice the initial stage where the migration rate decreases in the simulations with low viscosity ($\alpha_{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 - \dot{r}_P \tau_{memory} < r_0$). A similar behaviour can actually be seen for the 2D simulations, which we can now attribute to the gradual increase in the

$\mathcal{I}_{v,lib}/\mathcal{I}_{v,cross}$ ratio. For higher turbulent viscosities ($\alpha_{vis} \ge 10^{-4}$), $\mathcal{I}_{v,lib}/\mathcal{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 $\mathcal{I}_{v,lib}/\mathcal{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 $\alpha_{vis} = 10^{-4}$ compared to $\alpha_{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 ($\beta > 3/2$), we see that good agreement is obtained between the 1D and 2D simulations for $\alpha_{vis} \ge 10^{-4}$. For our lowest viscosity ($\alpha_{vis} = 10^{-5}$), our 1D model tends to underestimate $\mathcal{I}_{v,lib}/\mathcal{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 ($\tau_{vis}$) compares with the migration timescale ($\tau_{mig}$) and the libration timescale ($\tau_{lib}$) through the ratios $\tau_{vis}/\tau_{mig}$ and $\tau_{vis}/\tau_{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. $\dot{r}_P \ll 3\Omega_P x_s^2 / \{4\pi r_P\}$) so that Eq. (1) can be recast as (Paardekooper 2014)

$$ \Gamma_{c,dyn} = 2\pi \left( 1 - \frac{\mathcal{I}_{v,lib}}{\mathcal{I}_{v,cross}} \right) \Sigma_P r_P^2 x_s \Omega_P \dot{r}_P. \quad (14) $$

<footnote n="3">³ Recall that the torque formulae in Paardekooper et al. (2011) are for a smoothing length-to-scaleheight ratio $b/h = 0.4$, while our 2D simulations employ $b/h = 0.6$. The dependence of the individual torques on $b/h$ can be found in Paardekooper et al. (2010).</footnote>
Writing the total disc torque $\Gamma$ as $\Gamma_{\text{lin}} + \Gamma_{\text{c,static}} + \Gamma_{\text{c,dyn}}$ and as $\frac{1}{2} M_P r_P \Omega_P \dot{r}_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_P} \left(1 - \frac{\mathcal{I}_{\text{v,lib}}}{\mathcal{I}_{\text{v,cross}}}\right)}, \quad (15) $$

where $M_{\text{HS}} = 4\pi r_P x_s \Sigma_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 (\Gamma_{\text{lin}} + \Gamma_{\text{c,static}}), \quad (16) $$

with $f_{\text{red}} < 1$ being a reduction factor that reads

$$ f_{\text{red}} = \left[ 1 - \frac{M_{\text{HS}}}{M_P} \left( 1 - \frac{\mathcal{I}_{\text{v,lib}}}{\mathcal{I}_{\text{v,cross}}} \right) \right]^{-1}. \quad (17) $$

When $\Gamma_{\text{c,dyn}} \to 0$, $f_{\text{red}} \to 1$, whereas when $\Gamma_{\text{c,dyn}}$ is large, $f_{\text{red}} \ll 1$. Again, for a given disc model, $f_{\text{red}}$ depends on the planet's mass and orbital radius.

We evaluated the reduction factor $f_{\text{red}}$ for disc models similar to those used in global models of planet formation and evolution:

$$ \Sigma(r) = \Sigma_{5\text{au}} \times (r/5\text{au})^{-\beta}, \quad (18) $$

$$ T(r) = T_{1\text{au}} \times (r/1\text{au})^{-1/2}, \quad (19) $$

with $\Sigma_{5\text{au}} = 300 \text{ g cm}^{-2}$, $\beta \in [1/2, 1]$, and $T_{1\text{au}} = 280 \text{ 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, $\mathcal{I}_{\text{v,cross}}$ is obtained with Eq. (2) and $\mathcal{I}_{\text{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_0 < \dot{r}_P \cdot \tau_{\text{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 + \Gamma_{\text{c,dyn}}/(\Gamma_{\text{lin}} + \Gamma_{\text{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_{\text{red}}$ as a function of planet mass and orbital distance for four disc models with $\alpha_{\text{vis}} \in [10^{-5}, 10^{-4}]$ and $\beta \in [1/2, 1]$. Each panel shows two regimes, detailed below.

(i) $\tau_{\text{lib}}/2 > \tau_{\text{vis}}$ (below the solid white line): $f_{\text{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 $\mathcal{I}_{\text{v,lib}} \approx \mathcal{I}_{\text{v}}(r_P)$, and from Eq. (2) $\mathcal{I}_{\text{v,lib}}/\mathcal{I}_{\text{v,cross}} \approx 1 + (x_s/r_P)(3/2 - \beta)$. Eq. (17) then yields

$$ f_{\text{red}} \approx \frac{1}{1 + 4\pi\Sigma x_s^2 (\frac{3}{2} - \beta)/M_P}, \quad (20) $$

which is indeed independent of $\alpha_{\text{vis}}$ and $M_P$ (since $x_s \propto \sqrt{M_P}$) and decreases with radial distance for our radial profiles of density and temperature.

(ii) $\tau_{\text{lib}}/2 < \tau_{\text{vis}}$ (above the solid white line): here, a memory effect sets in for the librating vortensity and we see that $f_{\text{red}}$ now decreases with planet mass, but increases with alpha viscosity. We interpret this as a consequence of $\mathcal{I}_{\text{v,lib}}$ scaling as $(1 - \tau_{\text{memory}}/\tau_{\text{mig}})^{3/2-\beta}$ according to Eq. (6) with $\tau_{\text{memory}} \propto \tau_{\text{vis}}$ (Eq. (13)). Since $\tau_{\text{mig}} \propto 1/M_P$ and $\tau_{\text{vis}} \propto M_P^2/\alpha_{\text{vis}}$, increasing the planet mass or decreasing the alpha viscosity leads to a higher $\mathcal{I}_{\text{v,lib}}/\mathcal{I}_{\text{v,cross}}$ ratio, a stronger dynamical corotation torque, and, thus, a smaller reduction factor, $f_{\text{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 $\alpha_{\text{vis}}$ and lower $\beta$ (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 it 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_{\text{thermal}} = h^3 M_\star$. Another is the feedback mass, which is given by $M_{\text{feedback}} \approx 3.8 \times (Q/h)^{-5/13} \times 2/3 M_{\text{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 ($\ge 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, N_\phi) = (782, 2346)$, with grid cells spaced linearly or logarithmically along the radial direction. The full
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Maps showing reduction factor $f_{\text{red}}$, defined by Eq. (16), as a function of planet mass and orbital distance. The quantity $f_{\text{red}}$ 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 $f_{\text{red}} \to 1$ and strong where $f_{\text{red}} \ll 1$. Maps are obtained for four disc models that differ by the choice of the alpha turbulent viscosity ($\alpha_{\text{vis}}$) and the slope of the density profile ($-\beta$). Contour lines are plotted for $f_{\text{red}} \in [0.3, 0.5, 0.7, 0.9]$ for better visualisation. The solid white line marks the transition where $\tau_{\text{lib}}/2 = \tau_{\text{vis}}$ (see text). The dashed black line marks the feedback mass (see text). Regions where the planet mass exceeds the thermal mass (equivalent to $q_p/h^3(r_p) > 1$) are excluded as we expect substantial gap opening to occur in low-viscosity discs there.</figcaption></figure>

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_0 = 8, \beta = 1$ for $\alpha_{\text{vis}} = \{0, 10^{-7}, 10^{-6}\}$. The

simulations with $\alpha_{\text{vis}} = 10^{-7}$ and $\alpha_{\text{vis}} = 0$ show identical results that both deviate from those of the run with $\alpha_{\text{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 ($x_S \propto r_p^{1-f/2}$) 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_0 = 8$ runs for $\beta = \{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 $\dot{I}_{v,lib}$ and $\dot{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) \propto r^{-1}$ (the disc's flaring index thus cancels out) for $\alpha_{vis} \in \{10^{-5}, 10^{-4}\}$ and $\beta \in \{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 $\dot{I}_{v,lib}/\dot{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).

## 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 ($\dot{I}_{v,lib}$) and the orbit-crossing flow ($\dot{I}_{v,cross}$). Its inclusion in 1D simulations of planet formation and evolution necessitates expressions for $\dot{I}_{v,lib}$ and $\dot{I}_{v,cross}$, which is relatively straightforward for $\dot{I}_{v,cross}$ (Eq. (2)), but much less straightforward for $\dot{I}_{v,lib}$ since the latter depends on the planet's migration history. We show that the time evolution of $\dot{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 $\dot{I}_{v,lib}$ and $\dot{I}_{v,cross}$, which can be expressed as follows:

$$ \omega_{\text{lib}}(t) = \frac{\tau_{\text{vis}}\omega_0(r_p(t=0)) + t\omega(r_p(t))}{\tau_{\text{vis}} + t}, \qquad (21) $$

$$ \omega_{\text{cross}}(t) = \omega(r_p(t) - x_s(t)), \qquad (22) $$

where $\omega_0$ and $\omega$ denote the initial and instantaneous vortensity profiles, $\tau_{vis}$ is the viscous timescale across the planet's HS region (see Section 2.2), and $\dot{I}_{v,lib}(t) = \omega_{\text{lib}}^{-1}(t)$ and $\dot{I}_{v,cross}(t) = \omega_{\text{cross}}^{-1}(t)$. Eq. (22) is simply the same as Eq. (2). As for Eq. (21), the rationale behind the linear combination is as follows. When $t \ll \tau_{vis}$ (i.e. for very low turbulent viscosities), $\omega_{lib}(t) \approx \omega_0(r_p(t=0))$, the librating flow keeps its initial vortensity, as expected when vortensity advection dominates over vortensity diffusion. When $t \gg \tau_{vis}$, $\omega_{lib}(t) \approx \omega(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 $\dot{I}_{v,lib}/\dot{I}_{v,cross}$:

$$ \frac{\dot{I}_{v,lib}}{\dot{I}_{v,cross}} \approx 1 - \frac{x_s(t)}{r_p(t)} \frac{d \log \omega}{d \log r} \Big|_{r_p(t)} \quad \text{for } \tau_{\text{vis}} \ll t, \qquad (23) $$

$$ \approx \frac{\omega(r_p(t) - x_s(t))}{\omega_0(r_p(t=0))} \quad \text{for } \tau_{\text{vis}} \gg t. \qquad (24) $$

In the first limit (Eq. (23)), $\dot{I}_{v,lib}/\dot{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 $\dot{I}_{v,lib}$ and $\dot{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 $\beta < 3/2$, the alternative model better reproduces the initial increase in $\dot{I}_{v,lib}/\dot{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 $\beta > 3/2$ show quite remarkable agreement with the 2D simulations.

## 6.3. Comparison to previous work

A simple model for evaluating $\dot{I}_{v,lib}/\dot{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{\dot{I}_{v,lib}}{\dot{I}_{v,cross}} = \left(\frac{3}{2} - \beta\right) \frac{x_s^2}{6r_p\nu} j_p, \qquad (25) $$

which can be recast as

$$ \frac{\dot{I}_{v,lib}}{\dot{I}_{v,cross}} = 1 - \left(\frac{3}{2} - \beta\right) \frac{\tau_{\text{vis}}}{6\tau_{\text{mig}}}. \qquad (26) $$

When $\tau_{lib} \ll \tau_{vis} \ll |\tau_{mig}|$, the memory timescale in Eq. (13) is $\approx 0.4\tau_{vis} = 2\tau_{vis}/5$, and Eq. (6) yields

$$ \frac{\dot{I}_{v,lib}}{\dot{I}_{v,cross}} \approx 1 - \left(\frac{3}{2} - \beta\right) \frac{2\tau_{\text{vis}}}{5\tau_{\text{mig}}}, \qquad (27) $$

which is similar to the original model of Paardekooper (2014), except for the constant factor in front of $\tau_{vis}/\tau_{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 $\tau_{lib} \ll |\tau_{mig}|$, and using the expression in Eq. (26) for $\dot{I}_{v,lib}/\dot{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 $\dot{I}_{v,lib}/\dot{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 $\tau_{vis}/\tau_{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 $\dot{I}_{v,lib}/\dot{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 $\beta > 3/2$. This is not surprising since the above formulation for $\dot{I}_{v,lib}/\dot{I}_{v,cross}$ assumes a steady migration rate, which
<figure><img src="image_7.png" /><figcaption>**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.</figcaption></figure>

does not apply to the runaway migration regime. We briefly note that, due to a fast increase in the migration rate for $\alpha_{vis} = 10^{-5}$ and $\beta = 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) hydrodynamical 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 $\tau_{mig}$ is the characteristic timescale for the relative migration between the planet and disc. They predicted that super-Earths would have a bimodal
<figure><img src="image_8.png" /><figcaption>**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 ($\beta > 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.</figcaption></figure>

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
<figure><img src="image_9.png" /><figcaption>**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.</figcaption></figure>

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.

## 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 vortencies 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 vortencies 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 hydrodynamical 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 $\alpha_{vis} = 10^{-4}$, it is by a factor $\approx 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 ($\alpha_{vis} \lesssim 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.

## References

Alibert, Y., Mordasini, C., Benz, W., & Winisdoerffer, C. 2005, A&A, 434, 343

Bai, X.-N. 2013, ApJ, 772, 96

Bai, X.-N., & Stone, J. M. 2013, ApJ, 769, 76

Baruteau, C., & Masset, F. 2008, ApJ, 678, 483

Baruteau, C., & Masset, F. 2013, in Lecture Notes in Physics, eds. J. Souchay, S. Mathis, & T. Tokieda (Berlin: Springer Verlag), 861, 201

Baruteau, C., Bai, X., Mordasini, C., & Mollière, P. 2016, Space Sci. Rev., 205, 77

Benítez-Llambay, P., & Masset, F. S. 2016, ApJS, 223, 11

Benítez-Llambay, P., & Pessah, M. E. 2018, ApJ, 855, L28

Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112

Burn, R., & Mordasini, C. 2024, arXiv e-prints [arXiv:2410.00093]

Casoli, J., & Masset, F. S. 2009, ApJ, 703, 845

Chrenko, O., Chametla, R. O., Masset, F. S., Baruteau, C., & Brož, M. 2024, A&A, 690, A41

Crida, A., Griveaud, P., Lega, E., et al. 2022, in SF2A-2022: Proceedings of the Annual meeting of the French Society of Astronomy and Astrophysics, eds. J. Richard, A. Siebert, E. Lagadec, N. Lagarde, O. Venot, J. Malzac, J. B. Marquette, M. N'Diaye, & B. Briot, 315

de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529

Drążkowska, J., Bitsch, B., Lambrechts, M., et al. 2023, ASP Conf. Ser., 534, 717

Emsenhuber, A., Mordasini, C., Burn, R., et al. 2021, A&A, 656, A69

Emsenhuber, A., Burn, R., Weder, J., et al. 2023, A&A, 673, A78

Guilera, O. M., Benítez-Llambay, P., Miller Bertolami, M. M., & Pessah, M. E. 2023, ApJ, 95, 97

Hartmann, L., Calvet, N., Gullbring, E., & D'Alessio, P. 1998, ApJ, 495, 385

Kimmig, C. N., Dullemond, C. P., & Kley, W. 2020, A&A, 633, A4

Kimura, T., & Ikoma, M. 2022, Nat. Astron., 6, 1296

Lau, T. C. H., Birnstiel, T., Drążkowska, J., & Stammler, S. M. 2024, A&A, 688, A22

Lesur, G., Flock, M., Ercolano, B., et al. 2023, ASP Conf. Ser., 534, 465

Lin, D. N. C., & Papaloizou, J. C. B. 1993, in Protostars and Planets III, eds. E. H. Levy, & J. I. Lunine (Tucson: University of Arizona Press), 749

Liu, B., Lambrechts, M., Johansen, A., & Liu, F. 2019, A&A, 632, A7

Masset, F. 2000, A&AS, 141, 165

Masset, F. S. 2002, A&A, 387, 605

Masset, F. S. 2017, MNRAS, 472, 4204

Masset, F. S., & Papaloizou, J. C. B. 2003, ApJ, 588, 494

McNally, C. P., Nelson, R. P., Paardekooper, S.-J., Gressel, O., & Lyra, W. 2017, MNRAS, 472, 1565

McNally, C. P., Nelson, R. P., & Paardekooper, S.-J. 2018, MNRAS, 477, 4596

McNally, C. P., Nelson, R. P., Paardekooper, S.-J., & Benítez-Llambay, P. 2019, MNRAS, 484, 728

McNally, C. P., Nelson, R. P., Paardekooper, S.-J., Benítez-Llambay, P., & Gressel, O. 2020, MNRAS, 493, 4382

Mordasini, C., Alibert, Y., Klahr, H., & Henning, T. 2012, A&A, 547, A111

Müller, T. W. A., Kley, W., & Meru, F. 2012, A&A, 541, A123

Muto, T., & Inutsuka, S.-i. 2009, ApJ, 701, 18

Ndugu, N., Bitsch, B., Morbidelli, A., Crida, A., & Jurua, E. 2021, MNRAS, 501, 2017

Ogihara, M., Morbidelli, A., & Guillot, T. 2015, A&A, 584, L1

Ogihara, M., Kokubo, E., Suzuki, T. K., Morbidelli, A., & Crida, A. 2017, A&A, 608, A74

Ogihara, M., Kokubo, E., Suzuki, T. K., & Morbidelli, A. 2018, A&A, 615, A63

Paardekooper, S. J. 2014, MNRAS, 444, 2031

Paardekooper, S. J., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401, 1950

Paardekooper, S. J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293

Pierens, A. 2015, MNRAS, 454, 2003

Pinte, C., Teague, R., Flaherty, K., et al. 2023, ASP Conf. Ser., 534, 645

Rafikov, R. R. 2002, ApJ, 572, 566

Shakura, N. I., & Sunyaev, R. A. 1973, IAU Symp., 55, 155

Wafflard-Fernandez, G., & Baruteau, C. 2020, MNRAS, 493, 5892

Wafflard-Fernandez, G., & Lesur, G. 2023, A&A, 677, A70

Yun, H.-G., Kim, W.-T., Bae, J., & Han, C. 2022, ApJ, 938, 102

Zhu, Z., Stone, J. M., & Rafikov, R. R. 2012, ApJ, 758, L42

Ziampras, A., Nelson, R. P., & Paardekooper, S.-J. 2024a, MNRAS, 532, 351

Ziampras, A., Nelson, R. P., & Paardekooper, S.-J. 2024b, MNRAS, 528, 6130

Ziampras, A., Nelson, R. P., & Paardekooper, S.-J. 2025, arXiv e-prints [arXiv:2502.18564]
**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 = 0.05$. Since our fiducial disc model assumes $h(r) \propto r^{1/2}$, our simulations automatically tested the impact of having a smaller aspect ratio through migration (down to $h(0.5r_0) \sim 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 = 0.07$. We have increased the planet mass to maintain the same $q_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 $\Sigma_0$ to retain $Q_0 = 8$. The initial migration rate, which scales as $q_p/\{h_0Q_0\}$, is about twice as large as in our fiducial model for the same $Q_0$. These additional simulations used $\alpha_{vis} = 10^{-5}$ and $\beta = \{\frac{1}{2}, 1\}$. 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, $\bar{I}_{v,lib}/\bar{I}_{v,cross} - 1$ scales with $\tau_{vis}/\tau_{mig}$ (see Eqs. 26 and 27; these are valid for $\tau_{lib} \ll \tau_{vis} \ll |\tau_{mig}|$, and the latter inequality is marginally fulfilled in the simulations presented here). Since $\tau_{vis} \equiv q_p/\{\alpha h^3(r_p)\} \times \Omega^{-1} \propto r_p^{3/2-3f}$, $\tau_{vis}$ remains stationary for $f = 1/2$. When considering only the Lindblad torque for simplicity, $|\tau_{mig}| \propto r^{1/2+\beta}$ for $f = 1/2$: $|\tau_{mig}|$ decreases as the planet migrates inwards. This explains why $\bar{I}_{v,lib}/\bar{I}_{v,cross}$ also increases with time. Similarly, in the alternative model of Section 6.2, so long as $t \ll \tau_{vis}$, $\bar{I}_{v,lib}/\bar{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 \gg \tau_{vis}$. At the end of the simulations, $\tau_{vis}/t$ does not go below $\approx 0.5$. The final value for $\bar{I}_{v,lib}/\bar{I}_{v,cross}$ is larger in the $h_0 = 0.07$ simulations than in the $h_0 = 0.05$ runs because the planet has migrated farther inwards.

<figure><img src="image_10.png" /><figcaption>**Fig. A.1.** Direct comparison of 2D and 1D simulations results for $h_0 = 0.07$, $q_p = 2.7 \times 10^{-5}$ and $Q_0 = 8$. Rows are the same as in Fig. 2. Columns correspond to different approaches for calculating the inverse vortensities $\bar{I}_{v,lib}$ and $\bar{I}_{v,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.</figcaption></figure>
**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.

<figure><img src="image_11.png" /><figcaption>**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 $\alpha_{vis} = 10^{-3}$ are only shown for $\beta \ge 1.5$ to make plots less crowded.</figcaption></figure>