© The Authors 2023

1 Introduction

After more than 25 yr of observations, we have detected over 5000 exoplanets and have access to stunning observations of protoplanetary discs thanks to instruments including the Atacama Large Millimeter/sub-millimeter Array (ALMA). The exoplanet sample, together with our own Solar System, have revealed the existence of planets in our Galaxy of vastly different sizes, orbits, and bulk compositions, going from the smallest and densest terrestrial-type planets, to the super-Earths and mini-Neptunes with their moderate atmospheres, up to the giant Jupiter-like planets. In the meantime, high-angular-resolution images have recently shown the existence of detailed substructures in protoplanetary discs, such as rings, gaps, and spirals (e.g. DSHARP survey, Andrews et al. 2018; Huang et al. 2018, and MAPS programme, Öberg et al. 2021; Sierra et al. 2021), which challenge our understanding of the structure and evolution of such objects. Naturally, these two classes of astronomical objects are closely related to each other, so a crucial, bipartite question concerns how the protoplanetary disc environment sculpts the forming planetary system and simultaneously how forming planets shape the structure of the disc they are embedded within.

For one, a commonly proposed explanation for the substructures seen in DSHARP discs, such as rings and gaps in the millimetre emission together with corresponding deviations in Keplerian velocity of the gas, is that these features probe planet-disc interactions, thus revealing the ongoing process of planet formation (Pinilla et al. 2012; Zhang et al. 2018; Bae et al. 2018; Teague et al. 2018; Pinte et al. 2018); in at least one case, the well-known PDS-70 system, two forming protoplanets have been detected (Keppler et al. 2018; Wagner et al. 2018; Haffert et al. 2019). Even more so, many of these putative forming planets are found in a region of the orbital-period versus planetary mass planet that is currently unavailable to planet-detection methods used so far, due to observational biases (Lodato et al. 2019; Ndugu et al. 2019; Bae et al. 2018; Müller-Horn et al. 2022). The study of planet-disc interactions is therefore a useful tool to uncover exoplanetary systems from disc observations while they are forming. On the other hand, one of the main goals of planetary science is to construct a model to predict, at least in a statistical sense, what planetary system will emerge around a given star from a given disc (e.g. having given surface density, temperature, and thickness profiles, a given level of turbulent viscosity, etc; e.g. Ida & Lin 2008; Mordasini et al. 2009; Alibert et al. 2013; Alessi et al. 2017; Izidoro et al. 2017, 2021; Ndugu et al. 2018; Bitsch et al. 2019; Guilera et al. 2019; Emsenhuber et al. 2021). Clearly, such a model must take into account discplanet interactions, and their outcome can vary widely under different underlying disc environments.

Although they are unknown in our Solar System, the most common types of exoplanets in the galactic planetary census are the so-called super-Earths/mini-Neptunes (Fressin et al. 2013). These planets have masses of a few to a few tens of Earth masses and are observed very close to their host stars, often in compact multi-planetary configurations (Mayor et al. 2011; Winn & Fabrycky 2015; Weiss et al. 2018). Another relevant characteristic of these types of planets is that, when they appear in multi-planetary systems, their dynamical history is thought to have been shaped by mean motion resonant capture; most of these resonances have subsequently been broken after the disappearance of the disc (Izidoro et al. 2017, 2021; Pichierri & Morbidelli 2020; Goldberg et al. 2022), although a few remained stable in resonance and are still observed today (Trappist-1, Kepler-80, and Kepler-223, among others). Because of their abundance and their importance in planet formation theories, we concentrate on super-Earths in this work. Their planet-disc interactions typically fall in the so-called type-I migration regime, where the profile of the gaseous disc around them is only slightly perturbed from the underlying disc structure; the type-II regime instead pertains to planets that are massive enough to open a noticeable gap around their orbit, thus modifying the gas’ structure significantly. We concentrate here mainly on the type-I regime, although we also investigated the transition between type-I and type-II migration.

The study of planet-disc interactions in the type-I regime is the subject of a vast number of works, both analytical (e.g. Goldreich & Tremaine 1979, 1980; Artymowicz 1993; Ward 1997; Masset 2008) and numerical (e.g. Tanaka et al. 2002; Tanaka & Ward 2004; Kley & Crida 2008; Kley et al. 2009; Paardekooper et al. 2010, 2011; Jiménez & Masset 2017). Numerical studies resort to hydrodynamical simulations where a gaseous disc is fully resolved, to investigate the response of the planet to the perturbations driven onto the disc by the presence of the planet itself. This response usually manifests itself as a shrinking of the planetary orbit (a process called inward migration) and a damping of the orbit’s eccentricity and inclination (although outward migration and eccentricity excitation can also occur depending on the planet and disc parameters, Papaloizou et al. 2001; Kley & Dirksen 2006; Bitsch et al. 2013; Paardekooper & Mellema 2006; Bitsch & Kley 2010, 2011). Although they are extremely important tools, for example when one wants to compare disc models to real observations, such simulations are relatively costly and are for this reason not directly employable in population synthesis works, where instead N-body integrations with fictitious forces that mimic planet-disc interactions are preferred (e.g. Mordasini et al. 2009; Alibert et al. 2013; Izidoro et al. 2017, 2021; Bitsch et al. 2019; Ogihara & Hori 2020; Emsenhuber et al. 2021). One of the most widely used prescriptions for type-I interactions is that of Cresswell & Nelson (2008), who fitted dissipative evolution of super-Earth-type planets embedded in a gaseous disc to derive orbital element damping timescales that can be used as efficient recipes in N-body integrations. Cresswell & Nelson (2008) ran 3D simulations, using what are nowadays considered high-viscosity values. This may be an issue for low-viscosity (or thinner) discs, even for lower mass super-Earth-like planets, as even these may start carving a partial gap. Indeed, more recent non-ideal magnetohydrodynamical (MHD) simulations revealed that discs are less viscous than previously thought. While the question of the effect of different diffusive strengths has been thoroughly investigated for what concerns the planetary migration speed (i.e. the torque, e.g. Paardekooper et al. 2011, which is widely used in the literature in combination with the prescription of Cresswell & Nelson 2008 for eccentricity and inclination damping or, more recently, Jiménez & Masset 2017), this is not so in the case of eccentricity damping. Besides, as we mention earlier in this paper, it is known that heavier planets, which carve a significant gap, may even undergo eccentricity excitation due to planet–disc interaction (Papaloizou et al. 2001; Kley & Dirksen 2006; Bitsch et al. 2013; Duffell & Chiang 2015). The question naturally arises as of how to properly model the transition from classical (moderately high-viscosity, thick disc, and low-mass planet) type-I eccentricity damping and type-II evolution.

We note that for extremely low viscosities and moderately massive planets, vortices start to appear (Fung et al. 2014; McNally et al. 2019), making a clear description of planet-disc interactions elusive in these cases. In fact, the analytical formulas for the torque and eccentricity damping found in the literature are for the case of an axisymmetric disc with a smooth profile (Paardekooper et al. 2011; Jiménez & Masset 2017). Hydro-dynamical simulations seem to imply that, in the limit of an inviscid disc, planet-disc interactions for super-Earths lead to vastly different outcomes than in viscous discs, in that they can hinder capture into resonance and produce completely different systems (McNally et al. 2019). Moreover, analytical formulas that describe planet-disc interactions are typically only valid for eccentricities up to values comparable to the vertical aspect ratio of the disc (Tanaka & Ward 2004). Due to these limitations, we do not strive to derive general results in an arbitrarily large parameter space, and we instead focus on conditions that are normally encountered and useful to the community in a practical sense.

In particular, the aim of this paper is to revisit orbital eccentricity damping formulas for levels of viscosity that are more ‘modern’ compared to what was used in Cresswell & Nelson (2008), just as Paardekooper et al. (2011) included the effects of viscous diffusion in the expression of the Paardekooper et al. (2010) unsaturated torque. The general goal is to better understand the transition between viscous- and inviscid-disc type-I planet-disc interactions. First of all, even low-mass planets that would normally be considered in the type-I regime can, in low-viscosity environments, open significant gaps; this can be of interest in the context of the observational techniques mentioned above (Dong et al. 2017). However, most importantly, we are specifically interested in the process of resonant capture (McNally et al. 2019), since these are the mechanisms that are thought to sculpt the dynamics of the super-Earth population (Izidoro et al. 2017, 2021). These results will then be directly applicable to exoplanetary population synthesis models. Since we are not considering inclined orbits here, it is natural to consider a 2D setup. We also note that hydrodynamical investigations of super-Earths (such as the already mentioned McNally et al. 2019, but also others such as Ataiee & Kley 2021) are run for locally isothermal, 2D discs, since full 3D hydrodynamical simulations for resonant capture are too resource-intensive to be practical. In order to make a fair comparison with these works, and as a first step in our investigation, we also consider 2D, locally isothermal discs in this paper, while future work will be devoted to 3D effects, as well as the expression for inclination damping in partial-gap-opening scenarios. Besides, it is well known that resonant capture in the most common resonances (which are mean motion resonance of first order in the eccentricities) does not involve the inclinations, and all known planetary systems in confirmed resonances are observed to be extremely flat (Trappist-1, Kepler-80, Kepler-223, etc.). Moreover, 2D hydrodynamical simulations can be made to reproduce the main features of 3D simulations (Kley et al. 2009, 2012).

The rest of the paper is organised as follows. In Sect. 2, we describe our disc model, while in Sect. 3 we detail our hydrodynamical simulations and the methods used to describe partial gap opening and to calculate orbital element damping timescales, with specific attention to their implementation in N-body codes. In Sect. 4, we describe the results of our hydrodynamical simulations, while in Sect. 5 we discuss their implications. In particular, we present a simple recipe to model type-I eccentricity damping that takes into account the effects of partial-gap opening and that is consistent with hydrodynamical simulations. Finally, we conclude and summarise our results in Sect. 6. In addition, we give a simple theoretical argument to understand our main conclusions and supplementary explanations on our methods in the appendix.

2 Disc model

Our hydrodynamical experiments simulate a gaseous 2D disc around a M_* = M_⊙ star, extending from 0.35 to 3.3 AU. We assumed a power-law profile for the surface density $\text{Σ}\left(r\right)={\text{Σ}}_0{\left(r/r_0\right)}^{-{\alpha }_{\text{Σ}}}$ and a constant aspect ratio of H/r = h = h₀ across the disc (flaring index β_f = 0). We parameterised the turbulent viscosity ν_t of the disc using the following well-known alpha prescription (Shakura & Sunyaev 1973): ν_t = α_tH²Ω_K; we assumed a constant value for α_t. For such a disc, the accretion rate, Ṁ_gas,accr, is given by Ṁ_gas,accr = 3πνΣ_gas = 3πα_th²r²Ω_KΣ_gas. Assuming a constant accretion rate, we set α_∑ = 0.5. We also assumed that the locally isothermal approximation is valid; that is, we prescribed a temperature profile $T\left(r\right)=T_0{\left(r/r_0\right)}^{-{\beta }_T}$, where ß_T = 1 − 2ß_f = 1 by our aspect ratio prescription. The locally isothermal prescription assumes that cooling timescales are extremely short, and it is chosen here to provide a fair comparison with other works (McNally et al. 2019; Ataiee & Kley 2021). Although more sophisticated thermodynamical assumptions are beyond the scope of this work, we note that they may lead to different outcomes (e.g. in terms of gap opening, Miranda & Rafikov 2019, 2020; Zhang & Zhu 2020; note, however, that e-damping is similar in isothermal and fully radiative discs: Bitsch & Kley 2010) and will be the subject of future studies.

While the surface density and temperature profiles were fixed, h and α_t were left as free parameters. In our simulations, we used h ∈ {0.04, 0.05, 0.06} and α_t ∈ {3.16 × 10⁻⁵, 1. × 10⁻⁴, 3.16 × 10⁻⁴, 1. × 10⁻³}. Inside such a disc, we added a planet at a distance of a_pl = 1 AU using a sinusoidal mass taper to smoothly increase its mass from an initial value of 0 to a final mass of m_pl/M_* ∈ {1 × 10⁻⁵, 3 × 10⁻⁵, 6 × 10⁻⁵}, corresponding to typical masses of super-Earths, over the course of 50 orbits. The planetary masses are always below the corresponding thermal mass m_th = h³ M_*, so the disc-planet tidal perturbation does not drive local non-linear shocks and can be treated linearly (Lin & Papaloizou 1986). The planet’s eccentricity divided by the (fixed) aspect ratio was chosen as e_pl/h ∈ {0, 0.2, 0.4, 0.6, 0.8, 1}. We did not consider higher values for the eccentricity of the planet because the analytical formulas that we wish to compare our results to start breaking down. In any case, single planets in such mass ranges have their eccentricities damped by the disc, so extremely large eccentricities are not expected; instead, when multiple planets interact in a disc, for example by capturing in resonance via convergent migration, the expected capture eccentricities are of the order of h (Pichierri et al. 2018). Thus, considering e_pl up to a value of ≃h does not result in considerable drawbacks. Along the simulation, the planet’s orbit (semi-major axis and eccentricity) was kept fixed. We thus have a disc and planet setup with three free parameters, namely h, α, m_pl, for each of which we let e_pl/h vary. Table 1 lists all these parameters with the values chosen in our simulations. For each value of h, α_t, and m_pl, we run a simulation with e_pl/h spanning the reported values. The whole set is comprised of 216 high-resolution simulations.

We notice that the system is centred on the star and indirect forces should be considered. Recent works have shown that indirect terms must be carefully taken into account (Zhu & Baruteau 2016, Crida et al., in prep.). The recommendation of Crida et al. (in prep.) is to apply indirect forces to all the elements that feel a direct gravitational force. The planet feels the indirect force due to its own gravity as well as that of the disc. The disc feels indirect forces from the planet. Finally, we did not account for the indirect forces of the disc onto itself since we did not consider the disc’s self-gravity.

We implemented this setup in the fargOCA code (fargo with Colatitude Added; Lega et al. 2014)¹. The code is based on the fargo code (Masset 2000) extended to 3D. The fluid equations are solved using a second order upwind scheme with a time-explicit-implicit multi-step procedure. The code is parallelised using a hybrid combination of MPI and Kokkos (Carter Edwards et al. 2014; Trott et al. 2022). Code units are G = M_* = 1, and the unit of distance r₀ = 1 is arbitrary when expressed in AU. We use 1024 grid cells with arithmetic spacing in radius (corresponding to a δr ≃ 0.002) and 3000 cells in azimuth for the full (0, 2π; corresponding to a δϕ ≃ 0.002). Even for the smallest planetary masses, this ensures that we are resolving six cells in a half horseshoe width of such planets, which is needed in order to properly resolve the co-rotation torque (Paardekooper et al. 2011; Lega et al. 2014). We checked the convergence of our results by halving and doubling the resolution, seeing no discernible difference in the outcome. We used a smoothing length for the potential of the planet of r_sm = ϵH(r) with ϵ = 0.6, which better reproduces 3D effects in 2D simulations (Müller et al. 2012); finally, we used evanescent boundary conditions (de Val-Borro et al. 2006).

3 Methods

3.1 Disc density profile for partial-gap opening planets

A planet orbiting a star inside its protoplanetary disc will affect the disc structure in different ways depending on the system’s parameter. If the planet’s mass is low, its gravitational effect on the disc is small enough that the disc structure does not change significantly from the background unperturbed disc profile. This is the so-called type-I regime of planet-disc interactions. More massive planets will instead exert a torque onto the disc that can overcome the restoring viscous torque from inside the gas itself, and the planet will carve a gap around its orbit. This is the so-called type-II regime of gap-opening planets. For magnetorotational-instability- (MRI) type levels of viscosities (α_t ≳ 10⁻³), the type-I regime usually holds for planets up to a few tens of Earth masses. However, for sufficiently thin discs and at sufficiently low viscosities, (for example in the MRI-dead zone, e.g. Turner et al. 2014 for a review), even a small planet may start opening a partial gap.

We thus investigate the depth of the gap carved by simulated planets in a disc, to establish a threshold down to which type-I interactions can be considered valid. Crida et al. (2006) described the gap opened by a planet on a circular orbit by balancing gravity, viscous and pressure torques (the latter ones originating from the evacuation by pressure supported waves of gravitational torques). They define (arbitrarily) that a planet has opened a gap when the equilibrium disc surface density is 10% of the unperturbed disc surface density, Σ_min/Σ₀ = 0.1, and derive a condition for gap opening given by $$\frac{3}{4}\frac{H}{R_{\text{H}}}+\frac{50}{qℛ}\lesssim 1,$$(1)

where R_H is the Hill radius of the planet, q = m_pl/M_*, and $ℛ=r_{\text{pl}}^2{\text{Ω}}_{\text{K,pl}}/v_{\text{t}}$ is the Reynolds number.

Kanagawa et al. (2018) gave a prediction for the value of Σ_min/Σ₀, that is of the gap depth produced by a planet on a circular orbit, depending on the physical parameters of the disc and the planet. They showed that $${\text{Σ}}_{\min }/{\text{Σ}}_0\simeq \frac{1}{1+0.04K},$$(2)

where $$K:=q^2{h}^{-5}{\alpha }_t^{-1}$$(3)

is a dimensionless parameter. Gyeol Yun et al. (2019) slightly improved this result, replacing Eq. (2) with Σ_min/Σ₀ ≃ 1/(1 + 0.046 K). Finally, Bergez-Casalou et al. (2020) considered the effect of gas accretion in shaping the gap profile and gap opening mass for giant planets (typically q ≳ 10⁻⁴); since we are again interested in type-I interactions, we do not consider significant gas accretion in this work and instead keep our planetary masses fixed in time after the initial mass taper ramp-up.

It is interesting to note that Kanagawa et al. (2018) show that the transition between type-I and type-II migration happens at values of K of a few tens, which corresponds to a gap of the order of 0.5; their results seem, however, to depend on the level of viscosity, with the transition happening at larger values of K for lower viscosities. Instead, Crida et al. (2006) again considers a gap opened (and thus transition into type-II) when the gap is 10%. We can thus consider, for our purposes, that the transition from type-I to type-II migration happens when the gap depth is of a few 10⁻¹, and we won’t investigate the evolution below this threshold.

Since the eccentricities are small, we consider, for each setup shown in Table 1, the circular case as the nominal case to derive the gap profile. This allows us to compare our results with the results from the literature. We run all our e_pi = 0 simulations up to t_max = 3000 orbits of the planets, as was done in Kanagawa et al. (2018). In the lowest viscosity case, we integrated for an additional 1000 orbits to ensure the reaching of a steady-state surface density profile. Indeed, we check that over the last 50 orbits of the planet, the surface density does not change by more than ~0.1–0.5% (with the longest time needed for convergence given by the lowest-viscosity cases); as a further test of convergence, we integrated four setups spanning all the different α_t values for an additional 3000 orbits and checked that the relative difference in surface density over the last additional 3000 orbits is less than ≃2%. We then consider the surface density contrast, ${\text{Σ}}_{t_{\max }}\left(r\right)/{\text{Σ}}_0\left(r\right)$, where Σ_t is the (azimuthally averaged) surface density of the disc at time t, as a function of the radial distance, r. We mark the minimum value, ${\min }_r\left[\left({\text{Σ}}_{t_{\max }}/{\text{Σ}}_0\right)\left(r\right)\right]$, of the surface density contrast at a time of t_max and define it as the gap depth (or take an average of the two minima found slightly exterior and interior to the planet’s orbit), which we denote with Σ_min/Σ₀. Our results are presented in Sect. 4.2.

3.2 Orbital elements damping timescales

Along the simulation, at constant time intervals, the fargOCA code outputs two sets of forces: the (direct) force felt by the planet on its fixed orbit from the disc and the force felt by the star from the disc. The second force would not have any effect if the disc were axisymmetric, but because of the gas’ response to the presence of the planet, the induced asymmetry causes a net force felt by the star. Since in our simulation the frame of reference is astrocentric, we are not in an inertial reference frame. This means that the force felt by the star from the disc will result in an indirect (fictitious) response force felt by the planet in the astrocentric reference frame. Thus, we need to add this force to the one describing the direct planet–disc interaction. This yields a force, F_disc→pl, which describes the sum of direct and indirect planet-disc interactions, and thus the true force felt by the planet in an inertial reference frame.

From this force, we obtain the torque, $$\Gamma :=r_{\text{pl}}\times F_{\text{disc}\to \text{pl,}}$$(4)

and the power, $$P:=F_{\text{disc}\to \text{pl}}\cdot {\upsilon }_{\text{pl}}.$$(5)

We note that, since we are on the plane, the vector Γ = (0, 0, Γ_z)^⊤; so, we can concentrate on the scalar quantity Γ := Γ_z = ‖Γ‖. For a planet on an eccentric (planar) orbit, both the torque Γ and the power P are needed in order to determine the response of its orbital elements to the force F_disc→pl (we use orbit-averaged forces, where the average is done over 20 points along the planet’s eccentric orbit). This can be easily done as follows.

By definition, $\text{Γ:=}\frac{\text{d}{ℒ}_{\text{pl}}}{\text{d}t}={\dot{ℒ}}_{\text{pl}}$, that is, the rate of change of the (orbital) angular momentum, which is given by $${ℒ}_{\text{pl}}:=\text{m}\sqrt{\mu a\left(1-e^2\right)},$$(6)

where 𝔐 = (m_plM_*)/(M_* + m_pl) ≈ m_pl is the reduced mass of the planet, μ = 𝒢(M_* + m_pl) ≈ 𝒢M_* is the reduced gravitational parameter, and a and e are the semi-major axis and eccentricity of the planet’s orbit. At the same time, by definition, $P:=\frac{\text{d}{E}_{\text{pl}}}{\text{d}t}={\dot{E}}_{\text{pl}}$, that is, the rate of change of the (orbital) energy, which is given by $$E_{\text{pl}}:=-\frac{\mu \text{m}}{2a}.$$(7)

By the fundamental equation $\text{Γ}={\dot{ℒ}}_{\text{pl}}$, the angular momentum of the planet evolves, and we can introduce a migration timescale τ_m defined as $$\frac{{\dot{ℒ}}_{\text{pl}}}{{ℒ}_{\text{pl}}}=:-\frac{1}{{\tau }_{\text{m}}};$$(8)

this is, however, not the semi-major axis evolution timescale, and in the case of eccentric orbits the evolution of the eccentricity must also be taken into account. Indeed, we also define $$\frac{\dot{a}}{a}=:-\frac{1}{{\tau }_a},$$(9) $$\frac{\dot{e}}{e}=:-\frac{1}{{\tau }_e},$$(10)

in order to represent the evolution timescales of the semi-major axis and of the eccentricity of the planet. We now need to express these timescales in terms of the torque, power, angular momentum, and energy of the planet.

Taking the time derivative of the angular momentum and the energy, we obtain the following: $$P={\dot{E}}_{\text{pl}}=\frac{\mu \text{m}}{2a^2}\dot{a}=E_{\text{pl}}\left(-\frac{\dot{a}}{a}\right)\equiv \frac{E_{\text{pl}}}{{\tau }_a}$$(11) $$\begin{array}{l}\Gamma ={\dot{ℒ}}_{\text{pl}}=\text{m}\sqrt{\mu }\left(\frac{1}{2}{a}^{-1/2}\dot{a}{\left(1-e^2\right)}^{1/2}+\frac{1}{2}{a}^{1/2}{\left(1-e^2\right)}^{-1/2}\left(-2e\dot{e}\right)\right)\hfill \\ =ℒ\left(\frac{1}{2}\left(\frac{\dot{a}}{a}\right)-C\left(e\right)\left(\frac{\dot{e}}{e}\right)\right)\equiv ℒ\left(-\frac{1}{2{\tau }_a}+\frac{C\left(e\right)}{{\tau }_e}\right),\hfill \end{array}$$(12)

with $$C\left(e\right)=\frac{e^2}{1-e^2},$$(13)

C(e) ≈ e² for small eccentricities. Thus, we have $${\tau }_{\text{m}}=-\frac{{ℒ}_{\text{pl}}}{\text{Γ}},$$(14) $${\tau }_a=\frac{E_{\text{pl}}}{P},$$(15) $${\tau }_e=C\left(e\right){\left(-\frac{1}{{\tau }_{\text{m}}}+\frac{1}{2{\tau }_a}\right)}^{-1}.$$(16)

We also see that for a circular orbit τ_m = 2τ_a.

3.3 Planet-disc interactions in N-body codes

N-body codes implement type-I migration using the timescales τ_m and τ_e to define accelerations onto the planet given by (Papaloizou & Larwood 2000): $$a_{\text{m}}=-\frac{{\upsilon }_{\text{pl}}}{{\tau }_{\text{m}}},$$(17) $$a_e=-2\frac{\left({\upsilon }_{\text{pl}}\cdot r_{\text{pl}}\right)r_{\text{pl}}}{r_{\text{pl}}^2{\tau }_e},$$(18)

where r_pl and υ_pl are the planet’s position and velocity. The first equation describes a change in angular momentum, that is, a torque, so that the angular momentum evolves according to Eq. (8); the second equation represents a force that has zero torque since it is a radial force, so it does not contribute to Γ but only to the power, and implements, for small es,² an orbital damping of the eccentricity as described by Eq. (10). The semi-major axis evolution described by Eq. (9) thus results from a combination of torque and e-damping with a timescale given by $${\tau }_a={\left(\frac{1}{{\tau }_{\text{m}}/2}+\frac{C\left(e\right)}{{\tau }_e/2}\right)}^{-1},$$(19)

that is, the equivalent of Eq. (16).

3.3.1 Type-I forces

Various formulas exist in the literature to implement fictitious type-I forces in this framework. Cresswell & Nelson (2008) give explicit expressions for τ_m and τ_e (and τ_i, the damping of the inclinations for non-planar orbits, not considered in this paper) by fitting the orbital evolution of an m_pl = 10 M_⊕ planet placed on a variety of initial configurations. They used 3D simulations with a relatively low resolution (fewer than two cells per half horseshoe width of the planet, while in order to resolve the corotation torque one needs at least six cells, Paardekooper et al. 2011) and a fixed viscosity of α = 5 × 10⁻³. Defining the typical type-I damping timescale τ_wave (Tanaka & Ward 2004), $${\tau }_{\text{wave}}=\left(\frac{M_{\odot }}{m_{\text{pl}}}\right)\left(\frac{M_{\odot }}{{\text{Σ}}_{\text{gas,pl}}{a}_{\text{pl}}^2}\right)h^4{\text{Ω}}_{K,\text{pl}}^{-1},$$(20)

their best fit yielded (in the planar case) $${\tau }_{\text{m,CN2008}}=\frac{2{\tau }_{\text{wave}}}{2.7+1.1{\alpha }_{\text{Σ}}}{h}^{-2}\times P_e,$$(21)

where $$P_e=\frac{1+{\left(\widehat{e}/2.25\right)}^{1.2}+{\left(\widehat{e}/2.84\right)}^6}{1-{\left(\widehat{e}/2.02\right)}^4}$$(22)

is a reduction factor due to the planet’s eccentricity and ê = e/h, while $${\tau }_{e\text{,CN2008}}=\frac{{\tau }_{\text{wave}}}{0.780}\left[1-0.14{\left(\frac{e}{H/r}\right)}^2+0.06{\left(\frac{e}{H/r}\right)}^3\right].$$(23)

Paardekooper et al. (2010) focused instead on the explicit expression for the total torque experienced by a planet on a circular orbit, coupling linear estimates of the Lindblad torque Γ_L to a non-linear model of corotation torques. The latter torques are split into a barotropic (or vortensity-driven) component and an entropic (or thermal) component. Unlike the Lindblad torques, corotation torques are prone to saturation, so only the Lindblad torque would remain without any restoring process. Thus, Paardekooper et al. (2011) studied the effects of diffusion in restoring the corotation torque, ΓC, whereby in the limit of very strong diffusion, the linear corotation torque can be recovered. Their formulas are thus viscosity-dependent, and for this reason they are more cumbersome than those of Cresswell & Nelson (2008), so we do not re-write them here. They capture the behaviour of the torque observed in their 2D numerical experiments with an error of up to 20%. Jiménez & Masset (2017) re-analysed the torque experienced by low- and intermediate-mass planets using 3D simulations, deriving improved formulas valid for planets that do not significantly deplete their coorbital region. In addition to the vortensity and entropic components (scaling, respectively, with the gradient of vortensity and entropy), they consider a temperature component (scaling with the temperature gradient) and a so-called viscous coupling term. By striving for accuracy rather than simplicity and splitting the corotation torque into different components, the formulas of Jiménez & Masset (2017) are significantly more cumbersome, and we again do not reproduce them here. We note, however, that in our locally isothermal case, only the vortensity and temperature components appear.

An equivalent of the aforementioned studies on the detailed characterisation of the viscosity-dependent torque onto a planet but in the case of an eccentric planet (i.e. a formula for the power including diffusive effects) has not yet been carried out. Still, population synthesis models need to include both a and e damping (and inclination damping, not considered in this work), since planet-planet interaction can excite the planets’ eccentricities, and thus the circular approximation is not adequate in a practical sense. Many works (Izidoro et al. 2017, 2021; Emsenhuber et al. 2021) have implemented instead a combination of the formula of Paardekooper et al. (2011) for the torque and the expression of Cresswell & Nelson (2008) for the eccentricity damping, Eq. (23). Moreover, the torque itself has to be modified to include an e dependency (Bitsch & Kley 2010, 2011; Cossou et al. 2013; Pierens et al. 2013; Fendyke & Nelson 2014); this dependency introduces a reduction of the Lindblad torque similar to the one found by Cresswell & Nelson (2008), $${\text{Δ}}_{\text{L}}=P_e^{-1},$$(24)

and a reduction of the corotation torque given by $${\text{Δ}}_{\text{C}}=\exp \left(-e/e_{\text{f}}\right),$$(25)

where e_f = 0.5h + 0.01 is defined in Fendyke & Nelson (2014). The total torque is thus $${\text{Γ}}_{\text{tot}}={\text{Δ}}_{\text{L}}{\text{Γ}}_{\text{L}}+{\text{Δ}}_{\text{C}}{\text{Γ}}_{\text{C}}.$$(26)

Figure 1 contains an example of the comparison of orbital damping timescales, from one of our hydrodynamical simulations and from the formulas from the literature (CN2008 for Cresswell & Nelson 2008, P2011 for Paardekooper et al. 2010, and JM2017 for Jiménez & Masset 2017). The migration timescale, and thus the torque, give good agreements with these analytical predictions within their typical errors (of the order of 20%). The e-damping timescale is instead more efficient in the case shown here compared to the prediction from Cresswell & Nelson (2008). We compile our results for all of our numerical sample in Sect. 4.3.

3.3.2 Transition into type-II migration regime

Although we do not deal with the type-II migration regime in this paper, it is instructive to see how the transition from type-I to type-II migration is usually handled in N-body codes. Indeed, even though it is conceptually easy to think in terms of separate migration regimes, recent developments in our understanding of disc structure and evolution have shown that the processes that drive turbulent viscosities, such as the MRI, might be quenched in large portions of the midplanes where planets form; the remaining hydro-instability-driven viscosities would therefore be much lower than expected, so that even a low-mass planet may start opening a partial gap.

The hydrodynamical simulations of Kanagawa et al. (2018) allowed them to express the timescale of type-II migration τ_mig,II to the type-I migration timescales modulated by the gap depth Σ_min/Σ₀: $${\tau }_{\text{m,II}}={\left(\frac{{\text{Σ}}_{\min }}{{\text{Σ}}_0}\right)}^{-1}{\tau }_{\text{m,I}}.$$(27)

This can be used to transition from the two regimes in the case of partial-gap opening planets. After the gap is considered to have been fully opened (e.g. from Crida et al. 2006 condition (1)), the type-II migration regime (Lin & Papaloizou 1986; Robert et al. 2018) can be considered fully operational. We stress again that, so far, for the type-I regime only the migration timescale is modulated by the opening of a partial gap in population synthesis works, but not the eccentricity damping. Thus, in Sect. 4.3 we attempt to obtain an equivalent adaptation of Eq. (27) for the eccentricity damping for partial-gap opening planets.

Fig. 1 Comparison of orbital elements damping efficiencies 1/τmig (upper panel), 1/τe (middle panel), and 1/τa (lower panel) obtained from the output of our hydrodynamical simulations and from analytical formulas in the literature (CN2008 for Cresswell & Nelson 2008, P2011 for Paardekooper et al. 2010, and JM2017 for Jiménez & Masset 2017). The figures here are for the case of a mpl = 1 × 10−5M* planet with α = 3 × 10−4 and h = 0.05. τmig is obtained from the torque felt by the planet from the disc and typically shows good agreement between hydro and analytical predictions within the latter’s margin of uncertainty (20%). The analytical expression for τe comes from Cresswell & Nelson (2008) or from our fit presented in Sect. 4.3. τa is obtained from τmig and τe using Eq. (19).

4 Results

4.1 Emergence of vortices

Of our setups listed in Table 1, not all lead to a stable steady-state disc profile, as some of them showed the emergence of vortices. We chose not to include setups where vortices occurred in the final analysis presented in the next sections (we also tested that vortices appear in these setups even when a smoother mass taper is used, namely a ramp-up time longer by a factor of ten). This is because, even though vortices might dissipate over time (this typically happens over a few hundred to a few thousand orbits in our simulations), such cases cannot be fully captured by a simple analytical expression for type-I planet-disc interactions. In any case, we will see that vortices appear in our simulations when gaps are significantly deep for migration to be considered of the type-II regime, which is beyond the scope of this work. We note that whenever the e_pl = 0 simulation did not show the emergence of a vortex, neither did the e_pl ≠ 0 runs with the same value of h, α_t, and m_pl. This is advantageous both in theory and in practice; our results on the emergence of a vortex do not depend on the planetary eccentricity (across the values considered in this work), and each valid setup without vortices will contribute a data point for each value of the eccentricity. Five setups, (h = 0.04, α_t = 10⁻⁴, m_pl/M_* = 6 × 10⁻⁵), (h = 0.04, α_t = 3.16 × 10⁻⁵, m_pl/M_*, = 6 × 10⁻⁵), (h = 0.04, α_t = 10⁻⁴, m_pl/M_*, = 6 × 10⁻⁵), (h = 0.05, α_t = 3.16 × 10⁻⁵, m_pl/M_* = 3 × 10⁻⁵), and (h = 0.05, α_t = 3.16 × 10⁻⁵, m_pl/M_* = 6 × 10⁻⁵), showed the emergence of a vortex. These correspond, as a rule of thumb, to the cases of a lower α_t or h value, or the highest m_pl value, as expected in a qualitative sense. These simulations allow us to obtain a quantitative constraint on the system’s parameters that would lead to the emergence of a vortex. We give this constraint in the next section in terms of the depth of the gap carved by the planet. The case of (h = 0.05, α_t = 10⁻³, m_pl/M_* = 1 × 10⁻⁵), where a vortex does not appear, is shown as an example in Figs. A.1 and A.2.

4.2 Partial gap opening at low viscosity and/or thin discs

We calculated the gap depth, Σ_min/Σ₀, after t_max = 3000 orbits of the planet (or for an additional 1000 orbits for the lower viscosity cases) when a steady-state was reached, as described in Sect. 3.1. Figure 2 shows the outcome of our simulations compared to the predicted gap depths from Eq. (2) of Kanagawa et al. (2018). We note that they considered viscosities down to α_t = 10⁻³, while we go as low as α_t = 3.16 × 10⁻⁵. In order to make a fairer comparison, we run additional hydro-dynamical simulations with the same m_pl and h ranges from Table 1 for α_t = 3.16 × 10⁻³, in the e = 0 case alone, with the sole purpose of checking our gap-depth results against those of Kanagawa et al. (2018). We see that for α_t = 1−3 × 10⁻³ our experiments give a good match, as expected, to the results from Kanagawa et al. (2018), well within the spread observed in their simulations around their fit. However, for lower and lower viscosities, we observe a significant difference between the predicted value and the observed one. In general, Formula (2) from Kanagawa et al. (2018) overestimates how deep a gap a given planet will carve. In Fig. 3, we show an example with m_pl/m_* = 3 × 10⁻⁵, α_t = 3.16 × 10⁻⁴, and h = 0.05; the observed gap depth after 3000 orbits is ≃0.9, while the predicted value is ≃0.7. We note that after 3000 orbits of the planet, the surface density changes by less than 0.1% over 50 orbits.

Another feature seen in Fig. 2 is that the lower-viscosity cases are represented by fewer points. This is because in some of the simulations (namely those with more massive planet and/or thinner discs), a vortex has appeared that does not allow us to draw any definitive conclusion. Based on our numerical experiments, we can describe the limit in the parameter space after which one can expect to see a vortex in such simulations as a function of the gap depth. Figure 4 shows a diagram where we label the outcome of each simulation in green when a vortex has not appeared and red when it has appeared. When a vortex has not appeared, we report both the observed gap depth (filled circle) and the predicted gap depth (unfilled circle) from Kanagawa et al. (2018); when a vortex has appeared, we cannot use the simulations to observe a gap depth, and thus we only report the predicted gap depth from Kanagawa et al. (2018). We see that, in both cases, vortices are expected to appear when the gap depth Σ_min/Σ₀ ≲ 0.25.

Fig. 2 Predicted gap depth, Σmin/Σ0, from Kanagawa et al. (2018) versus value obtained from the simulations. Different colours represent different levels of viscosity as shown in the legend. The lower the viscosity, the bigger the difference in gap depth. However, down to the viscosities considered by Kanagawa et al. (2018), we see good agreement, as expected (the data for αt = 3.16 × 10−3 are shown with smaller symbols as they are only used for the analysis of the gap depth). Different shapes are used to represent different aspect ratios and planetary masses; squares, downward-pointing triangles, and upward-pointing triangles represent aspect ratios of 0.04, 0.05, and 0.06, respectively; empty, filled and crossed symbols represent mpl/M* = 10−5, 3 × 10−5, and 6 × 10−5, respectively.

Fig. 3 Surface density profile, Σt/Σ0, after ${\overline{t}}_{\max }=2950$ orbits and tmax = 3000 orbits of a mpl/M* = 3×10−5 planet in a disc with αt = 3.16 × 10−4 and h = 0.05. We note that after the last 50 orbits, the surface density has reached a steady state, with a change of less than 0.1%. The black horizontal dashed line represents the prediction of the gap depth from Kanagawa et al. (2018). The two dot-dashed horizontal lines are used here to measure the observed depth of the gap carved by the planet by taking the average of their values.

Fig. 4 Emergence of vortices in our simulations as a function of the gap depth. When no vortex was observed, we report both the observed gap depth at the end of our simulations, as well as the predicted one from Kanagawa et al. (2018), with a green circle (filled and unfilled, respectively). A red ⊗ symbolises those simulations where a vortex appeared, in which case only the predicted value from Kanagawa et al. (2018) is plotted. A dashed vertical line marks the approximate location of the transition between no vortex and vortex.

4.3 Eccentricity damping efficiency for partial-gap opening planets

Following Sect. 3.2, we extract the eccentricity damping efficiency of 1/τ_e from our hydrodynamical simulations for all setups where no vortex has emerged. Figure 5 shows the eccentricity damping efficiency, 1/τ_e, normalised by the expected efficiency, 1/τ_e,CN2008, from Cresswell & Nelson (2008, Eq. (23)), which is the one that has been used extensively in the literature so far (e.g. Izidoro et al. 2017, 2021; Emsenhuber et al. 2021). We show the damping efficiency as a function of the gap depth carved by the planet (see previous sections), both using the prediction of Kanagawa et al. (2018, panel a) and the gap obtained from each hydrodynamical simulation in the circular case (panel b). While panel a yields a very noisy plot, it is clear in panel b that there is a strong trend when one uses the actual gap depth.

The overall trend is, however, clear in both plots; the eccentricity damping is less efficient for deeper gaps (left side of the plots), down to a factor of ~1/4 less efficient at gap depths of ≃0.25, which is close to the transition from type-I to type-II regimes and where we start obtaining vortices in our setups (Fig. 4). This suggests that, just as the transition from type-I to type-II migration speed is modulated by the gap depth (Eq. (27)), so should the eccentricity damping. We now look for a more quantitative expression of this fact, and since the data are much cleaner in panel b, we use the observed gap depth rather than the predicted one. We deal with the question of how to predict a gap depth in an N-body setting without resorting to computationally expensive hydrodynamical simulations later in the paper.

For low eccentricities, 0 < e/h ≲ 0.4, one can fit 1/τ_e versus the observed gap depth very well with a straight line over the full gap depth range considered here (Σ_min/Σ₀ ≃ 0.3 to 1). In the limit of no gap (Σ_min/Σ₀ ≃ 1), which should correspond to the setup of Cresswell & Nelson (2008), the observed eccentricity damping efficiency is slightly higher than the prediction 1/τ_e,CN2008 (or in other words, τ_e is slightly smaller than τ_e,CN2008). We note that we used a 2D setup with very high resolution, while Cresswell & Nelson (2008) used a 3D setup with lower resolution. Moreover, the hydrodynamical codes used are different, and Cresswell & Nelson (2008) fitted e-damping timescales to evolving planets, while we keep the planet on a fixed orbit and extract τ_e from the forces felt by the planet from the disc. We also ran cases with the same setup as Cresswell & Nelson (2008) and found that we still obtain a slightly more efficient e-damping. The difference, however, is smaller than 20%, which is the accuracy of the Paardekooper et al. (2011) torque formula anyway. For higher eccentricities up to e/h ≃ 1, the points again follow a straight line for gap depths ≃0.3 up to ≃0.8, after which e-damping becomes super-linear and significantly more efficient than the prediction of Cresswell & Nelson (2008). We checked this result with the original fargo code of Masset (2000), yielding very similar results to ours for the orbital damping timescales, with a difference of only 1%. One explanation for this effect is that for shallower gaps, and thus also thinner gaps, and sufficiently high e, the planet’s excursions around r_pl = a due to the eccentricity of the orbit start having a significant effect. Since the gap around a planet starts to be carved at r_pl ± 2/3H where the Lindblad torques accumulate, we see that the planet’s orbit starts interacting more with the edge of the gap at e/h ≃ 1. Similar effects were observed in Bitsch & Kley (2010) and Fendyke & Nelson (2014).

We thus consider a fit to the data separating the horizontal axis into gap depths deeper and shallower than 80%, and we perform a double fit over the two portions that join continuously at gap depths of 80% (Fig. 6). We find that the simple piecewise linear fit gives a very good approximation to the data: $$\begin{array}{l}\frac{1}{{\tau }_e}=\frac{1}{{\tau }_{e\text{,CN2008}}}\times \{\begin{array}{cc}{c}_1-m_1\left(1-\frac{{\text{Σ}}_{\min }}{{\text{Σ}}_0}\right)& \text{if}\frac{{\text{Σ}}_{\min }}{{\text{Σ}}_0}<0.8,\\ c_2-m_2\left(1-\frac{{\text{Σ}}_{\min }}{{\text{Σ}}_0}\right)& \text{if}\frac{{\text{Σ}}_{\min }}{{\text{Σ}}_0}>0.8\end{array}\hfill \\ m_1=1.263,c_1=1.169,m_2=\tilde{m},c_2=0.916+0.2\tilde{m},\hfill \end{array}$$(28)

where $\tilde{m}=\max \left\{m_1,m_1+3\left(e/h-0.3\right)\right\}$. In particular, the multiplicative factor would in general depend on e/h, but its expression is rather simple. In Fig. 6, we show a comparison between the fit and the data. The typical error given by the fit is of the order of 4%, and it is overall less than 20%, which is the accuracy of torque formulas from the literature. When Σ_min/Σ₀ = 1 (no measured gap), we recover the results of Cresswell & Nelson (2008) within errors equivalent to the error to the torque typically used in the literature (Paardekooper et al. 2011), while for deeper and deeper gaps Eq. (28) yields a simple model for the reduction of e-damping efficiency due to partial gap opening.

5 Discussion

5.1 Practical estimation of the gap profile

The results of Sect. 4.2 on the shape and depth of the gaps carved by planets indicate that a more careful analysis is needed when going to lower disc viscosities. In addition, for our own purposes of fitting e-damping timescales for use in N-body codes, the two panels from Fig. 5 show that a better fit for the gap depths at lower viscosities is needed in order to make any practical use of our Eq. (28). We recall that, down to α_t viscosities of about 10⁻³, the prediction of Kanagawa et al. (2018) also shows relatively good agreement with our numerical experiments (with a perfect fit at α_t = 10⁻³, see Fig. 2); moreover, Eq. (2) has the unquestionable advantage of being very easy to handle in practice for N-body codes and population synthesis studies. Thus, we take a practical approach and build upon the ideas of Kanagawa et al. (2018) to obtain a similar fit to our data. In particular, we rewrite Formula (2) from Kanagawa et al. (2018) as $$K=25\left(\frac{1}{{\text{Σ}}_{\min }/{\text{Σ}}_0}-1\right),$$(29)

and we fit a function of the form $\tilde{K}\left(q,{\alpha }_{\text{t}},h\right)=C{q}^{p_q}{\alpha }_t^{p_{{\alpha }_t}}{h}^{ph}$ using our observed values for the right hand side of Eq. (29). We used Mathematica’s NonlinearModelFit function and obtained the following fit: $$\tilde{K}=3.93q^{2.3}{h}^{-6.14}{\alpha }_t^{-0.66}.$$(30)

We show, in Fig. 7, an analogue of Fig. 2, but using $\tilde{K}$ from Eq. (30) instead of K from Eq. (3), showing much better agreement. This suggests a modification to the gap depth formula, namely $${\text{Σ}}_{\min }/{\text{Σ}}_0\simeq \frac{1}{1+0.04\tilde{K}}.$$(31)

We should stress that Eq. (30) is only a fit to the data and does not emerge from a careful study of the gas dynamics at low-to-moderate viscosities due to the presence of a planet. However, its sole purpose is to be combined with Eq. (28) in order to give a practical recipe for e-damping without resorting to hydrodynamical simulations (see next section). Moreover, we note that 2D and 3D simulations can be different regarding gap opening (see appendix of Bitsch et al. 2018), which was not considered here. A more detailed analysis on the matter will be the subject of future studies, as well as a study in the full 3D case. We note that Fung et al. (2014) and Bergez-Casalou et al. (2020) already showed that the gap opening at low viscosities does not follow the predictions from higher viscosity simulations α_t ≳ 10⁻³. Also, one should always keep in mind that these studies have been carried out using different hydrodynamical codes.

Fig. 5 Eccentricity damping efficiency versus gap-depth for all setups where no vortex emerged. In both panels, the e-damping efficiency on the vertical axis is normalised by the expected value from Cresswell & Nelson (2008). A 20% error around the expected value in the limit of no gap (to the right in the plots) is shown by two dashed horizontal grey lines. Values for the e-damping efficiency are shown for different eccentricities, e/h ∈ {0.2, 0.4, 0.6, 0.8, 1), by points of different colours joined together by opaque lines. Panel a shows on the horizontal axis the gap depth according to the prediction of Kanagawa et al. (2018), which results in an extremely noisy scatter plot. Panel b shows, on the horizontal axis, the gap depth observed from the actual simulations as plotted in Fig. 2, which shows a much cleaner dependence. In both panels, the overall trend clearly shows a decrease in e-damping efficiency (i.e. longer e-damping timescales) for deeper and deeper gaps, down to a factor of ~l/4 less efficient eccentricity damping at the transition from type-I to type-II regimes (gap depths of ≃0.25) as compared to the limit of no gap.

Fig. 6 Double linear fit of e-damping efficiency to the data from panel b of Fig. 5 according to Eq. (28). The fit to the data is split into a linear fit over gap depths ∈ [0.3, 0.8], and a separate linear fit over gap depths ∈ [0.8,1], which continuously joins with the previous fit. The slopes for both fits are given in the legend in the top left corner.

Fig. 7 Similar to Fig. 2, but using Eq. (31) as an approximation for the gap depth, showing better agreement with the data. As in Fig. 2, different colours represent different levels of viscosity as shown in the legend, and different symbols represent different aspect ratios and planetary masses: squares, downward-pointing triangles, and upward-pointing triangles for h = 0.04, 0.05, and 0.06 respectively; empty, filled, and crossed symbols for mpl/M* = 10−5, 3 × 10−5, and 6 × 10−5, respectively.

Fig. 8 Similar to Fig. 6, but using Eq. (31) as an approximation for the gap depth on the horizontal axis. The fit to the e-damping efficiency obtained from hydrodynamical simulations remains good, within the typical variations of analytical formulas for type-I migration.

5.2 Implementation of eccentricity damping to N-body codes for partial gap-opening planets

Using the approximate gap depth of Eq. (31), we obtain, together with Eq. (28), a practical recipe to implement e-damping in low viscosity discs for partial-gap-opening planets. The resulting comparison to the data, using Eq. (31) instead of the observed gap depth as in Fig. 6, is shown in Fig. 8. The typical error of this fit is still of the order of 4%, and less than the ≃20% of the torque formalisms. An example of comparison between the output of our hydrodynamical simulations and the fit is already shown in Fig. 1.

Concerning the left hand side of Fig. 8, we recall that according to Kanagawa et al. (2018) the transition between type-I and type-II migration regimes should occur when the gap depth is about 50% (see their Fig. 4); instead, Crida et al. (2006) considered (arbitrarily) that this transition should happen when the bottom density is ≲10% of the unperturbed surface density. Such deep gaps can be achieved for massive planets, but it is unclear whether a low-mass planet of a few (or a few tens) of Earth masses in a low-viscosity disc would be able to carve such gaps without the emergence of a vortex. In any case, we see in Sect. 3.1 that a gap of the order of a few times 10⁻¹ is the threshold down to which one can consider the type-I regime valid. This means that we are able to probe reasonably deep gaps, even in the lowest viscosity cases considered here, down to the type-I to type-II transition, so the results presented in this paper can be safely used in N-body implementations where a real disc is not evolved. Intermediate gap depths seem to be the easiest to grasp: The e-damping efficiency is simply linear with the gap depth, with very little dependence on the eccentricity (see also Appendix B). The e-independent fit reported in the first line of Eq. (28) captures the main effect of reduced e-damping efficiency for intermediate gaps with respect to the commonly used Eq. (23) from Cresswell & Nelson (2008); it is thus very similar to the expression for the change in migration timescale for partial-gap-opening planets, Eq. (27), which is the approach that has been typically used in the literature. Concerning instead the right hand side of Fig. 8, the deviation from the e-independent linear fit at high eccentricities, as mentioned in Sect. 4.3, seems to be caused by the interaction with the edge of the gap. The latter is known to show differences in 2D versus 3D hydro-dynamical simulations. Moreover, there are known 2D and 3D differences for the change of eccentricity (Bitsch & Kley 2010). Thus, some differences might be expected when going into the 3D case. This will be the subject of future work.

5.3 Consequences for mean motion resonant locking via convergent migration

There is strong theoretical and observational evidence that the dynamical histories of super-Earth/mini-Neptune systems have been shaped by a phase of type-I-migration-driven capture in mean motion resonance (Izidoro et al. 2017, 2021). Hydrodynamical simulations from McNally et al. (2019) suggest instead that, in the limit of an inviscid disc, resonant capture might be avoided by planets that would otherwise build a resonant chain in a viscous disc. Their setup includes five planets of increasing mass with increasing separation from the central star inside a 2D disc with a prescribed surface density and temperature profile, in the case of a viscous (α_t of the order of 10⁻³) and inviscid disc. While in the viscous case the stable systems were always found to be locked in mean motion resonances, the inviscid case did not form complete chains despite the final system being closely packed. McNally et al. (2019) did not give a definitive explanation for the difference in outcome for viscous and inviscid discs. They mention that significant perturbations to the disc structure can occur in the inviscid case (such as vortices), which may be the reason why some planets escape mean motion resonances and undergo chaotic close encounters. It is, however, not easy to draw definitive conclusions on exactly what this outcome is due to and whether it is representative of the bulk of the super-Earth/mini-Neptune population. Indeed, it is hard to determine a precise dynamical history from the evolution of a five-planet system; also, a setup with increasing planetary mass is in effect built to favour instabilities, since resonances tend to be less stable when the outer planet is more massive than the inner one (Morbidelli et al. 2008; Deck & Batygin 2015). A completely inviscid disc might not be a fully realistic environment, as even in the MRI dead zone a residual of turbulent viscosities can arise by purely hydrodynamical instabilities, with α_t of the order of 10⁻⁴ (e.g. Pfeil & Klahr 2021; Flock et al. 2020), and observational constraints determine α_t in discs to range between 10⁻⁴ to a few 10⁻² (Rafikov 2017; Dullemond et al. 2018; Flaherty et al. 2017, 2018). Finally, we do observe a few systems in resonance, and their observed architecture has been used to determine, a posteriori, the level of viscosity; in the case of Kepler-223, this indirect approach has yielded an estimate on α_t of the order of 10⁻³ (Hühn et al. 2021).

In any case, in light of the results of McNally et al. (2019), one should expect a transition from efficient to inefficient resonant chain assembly as planet-disc interactions become more and more significant even in the type-I regime, which is due to lower viscosity discs. Our results from Fig. 5 might give an initial hint as to why this is the case, even when vortices are not formed by these planets. Indeed, we found that the eccentricity damping provided by the disc can be significantly less effective for deeper and deeper gaps. During resonant capture, the orbital eccentricities are pumped by the planet-planet interaction until an equilibrium configuration is obtained between resonant pumping and disc-driven e-damping (Pichierri et al. 2018). Less efficient e-damping thus means higher orbital eccentricities, which lead to a more excited state more prone to instabilities. Finally, the emergence of vortices may indeed play a role, as McNally et al. (2019) suggest. Clearly, more work is needed in order to obtain a more detailed picture of type-I-like interactions for low-viscosity discs in the context of shaping the exoplanet population and for applications to population synthesis models.

We conclude remarking the work of Ataiee & Kley (2021), which performed 2D, locally isothermal simulations of resonant trapping for low- (α_t = 5 × 10⁻⁵) viscosity discs in a variety of configurations, implementing two different types of disc inner edge to stop the inner planet’s migration. Since they also considered increasing planetary masses for further-out planets, some simulations failed to build chains, remarking that the more massive the outer planet is, the earlier the system becomes unstable, which is in line with the aforementioned N-body and theoretical investigations on the subject. Still, their simulations with two-to-three planets were able in many cases to form resonant chains that remained stable also after the removal of the gas. Indeed, resonant systems of two-to-three planets are known to be relatively stable for low-mass planets and low eccentricities (Pichierri et al. 2018; Pichierri & Morbidelli 2020). For a higher number of planets or planetary cores, a less efficient e-damping in low-viscosity (but not inviscid) discs allows for chains to form during the disc phase, which, because of the increased eccentricities with respect to classical type-I migration e-damping, become more susceptible to instabilities after the dissipation of the disc. This would have consequences in population synthesis works, such as the so-called breaking the chains scenario (Izidoro et al. 2017, 2021), allowing instabilities to also develop for planets of lower masses than the ones needed in these works to efficiently break the resonant chains.

6 Conclusions

In this paper, we revisit the problem of orbital eccentricity damping timescales for planets that are in the type-I migration regime but open partial-to-significant gaps, for example due to the low viscosity in their surrounding disc. These planets, typically in the super-Earth/mini-Neptune mass range, are therefore not expected to undergo type-II migration, but their interaction with thin, low-viscosity discs can be viewed as pertaining to a transition regime between type-I and type-II migration. In general, the shape and depth of the gap carved by the planet and the resulting change in migration timescales (i.e. the torque felt by the planet) in such intermediate cases (up to and including classical type-II interactions) have been the subject of many theoretical and numerical works (e.g. Crida et al. 2006; Paardekooper et al. 2011; Jiménez & Masset 2017; Kanagawa et al. 2018); these theories have been implemented in multiple population synthesis works (e.g. Izidoro et al. 2017, 2021; Ndugu et al. 2018; Bitsch et al. 2019; Ogihara & Hori 2020; Emsenhuber et al. 2021) to reproduce the observed characteristics of exoplanetary systems. In particular, the transition in migration regimes from type-I and type-II has been described as a simple reduction to the type-I torque proportional to the gap-depth carved by the planet (Kanagawa et al. 2018). However, a similar study for the change in eccentricity damping has not been performed. Instead, the aforementioned population synthesis works have used a fixed e-damping efficiency following the work of Cresswell & Nelson (2008), which fitted the evolution of Super-Earths embedded in a disc of relatively high viscosity (α_t = 5 × 10⁻³) to obtain orbital damping timescales. These high α_t values might not, however, be the norm in planet forming regions, meaning that even such low-mass planets as the ones considered in Cresswell & Nelson (2008) would open partial gaps at lower viscosities (α_t ≲ 10⁻⁴). Besides, it is known that the orbital eccentricities can even be pumped by planet-disc interactions for planets that have carved a deep-enough gap (Papaloizou et al. 2001; Kley & Dirksen 2006; Bitsch et al. 2013). We thus expect there to be a transition between ‘pure’ type-I eccentricity damping as in Cresswell & Nelson (2008) and less efficient e-damping up to even e-pumping for planets that carve intermediate to deep gaps.

We thus explored this transition making use of high-resolution, 2D, locally isothermal hydrodynamical simulations of super-Earth-type planets of varying masses (m_pl/M_* = 1 to 6 × 10⁻⁵) and varying orbital eccentricities (e/h = 0 to 1, where h = H/r is the disc’s aspect ratio), embedded in discs of varying viscosities (α_t = 10⁻³ down to 3.16 × 10⁻⁵) and aspect ratios (h = 0.04 to 0.06). The planet is kept on a fixed orbit and the system is evolved for thousands of the planet’s orbital periods until a steady-state is achieved, all the while recording the planet–disc interaction forces that drive migration and e-damping. The e = 0 case is used in order to measure gap depths of ${\text{Σ}}_{t_{\max }}/{\text{Σ}}_0$, where ${\text{Σ}}_{t_{\max }}$ is the radial surface density profile at the end of the simulation and Σ₀ is the initial, unperturbed surface density (a power law). We compared our observed gap depths with the already existing results from Kanagawa et al. (2018), showing good agreement at the higher viscosities considered by Kanagawa et al. (2018, down to α_t ≃ 10⁻³), but less so for the lower viscosities; in particular, their formula gave an overestimation of the gap depths observed in our own simulations.

In Sect. 5.1, we give a fit to the gap depth similar to that of Kanagawa et al. (2018), which shows better agreement with our low-viscosity simulation, as a function of the disc’s aspect ratio and viscosity and the planet’s mass. Then, we recorded the observed e-damping efficiency for the various planetary eccentricities, as a function of the observed gap depth. For deep gaps, a clear linear trend is observed, with e-damping efficiency scaling with the gap depth and being as low as 1/4 of the expected value from the prescription of Cresswell & Nelson (2008) at the transition between type-I and type-II migration. The result is thus similar to the case of diminished migration efficiency for partial-gap opening planets proposed in Kanagawa et al. (2018) and used extensively in the aforementioned population synthesis works. For shallow gaps, the same simple linear trend remains for low eccentricities (e/h ≲ 0.2), while we observe an increase in eccentricity damping efficiency for higher eccentricities. This can be understood as due to the interaction of the eccentric planet with strongly non-axisymmetric structures in the disc, in particular with the gap carved by the planet itself.

A simple e-dependent fit is found for all cases in Eq. (28), which yields typical errors smaller than or comparable to the accuracy of widely used, analytical type-I migration prescriptions (Paardekooper et al. 2011). Together with the gap-depth fit of Eq. (31), this allows us to provide a simple prescription for eccentricity-damping in the type-I regime for partial-gap opening planets that is consistent with high-resolution 2D hydrodynamical simulations. 3D simulations as well as more sophisticated thermodynamics to relax the locally isothermal assumption will be the subject of future work in order to obtain a full picture for eccentricity and inclination damping efficiency for partial-gap-opening planets.

Acknowledgements

G.P. and B.B. thank the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. We acknowledge HPC resources from “Mesocentre SIGAMM” hosted by Observatoire de la Côte d’Azur. L.E. thanks Alain Miniussi for maintenance and re-factorisation of the code FARGOCA. G.P. thanks Frédéric Masset for helpful discussions. The authors wish to thank the anonymous referee for helpful comments and suggestions.

Appendix A Disc surface density and torque response

We present snapshots from a typical evolution of one of our discs, with α_t = 10⁻³ and h = 0.05, interacting with an m_pl/M_* = 1 × 10⁻⁵ planet. Figure A.1 depicts the surface density contrast, showing the emergence of the wake typical of a type-I response of the disc due to the presence of the planet; Figure A.2 shows the stabilisation of the torque felt by the planet from the disc.

Fig. A.1 Surface density evolution in a setup with αt = 10−3 and h = 0.05, over 3000 orbits of a mpl/M* = 1 × 10−5 planet on a (fixed) circular orbit at 1 AU. We show the surface density contrast (∑t(r, θ) − 〈∑t(r)〉θ)/〈∑t(r)〉θ on the (θ, r) plane, where 〈∑t(r)〉θ is the azimuthally averaged surface density at the distance r.

Fig. A.2 Evolution of the torque felt by the planet in the setup presented in Figure A.1 as a function of time. The vertical dashed line represents timestamps for which the surface density is shown in Figure A.1. The torque stabilises quickly in this case, as the system reaches a long-lived steady state.

Fig. A.3 Linear fit of the e-damping efficiency to the data from panel b of Figure 5, that is, using the gap depths observed from the simulations on the horizontal axis. In panel a we fit the full data, both for the separate values of e/h (coloured dashed lines) and for the full data set combining all eccentricities (black dashed line). In panel b we fit the data in the same way, but considering gap depths of up to 80% (indicated by the dashed red vertical line). In both cases we report the observed slopes of the lines in the top left legend.

Appendix B Estimate of eccentricity-damping and pumping resonances for partial gaps

For an eccentric planet, the Fourier decomposition of the gravitational planet–disc potential includes terms that are proportional to the eccentricity and therefore drive its evolution. Typically, for a small planet that does not significantly perturb the disc profile, the eccentricity is both damped and pumped by these individual so-called first-order Lindblad resonances (although their combination usually results in damping of the eccentricity, see below). The main resonances that damp e are the coorbital first-order Lindblad resonances with |1 − m| = 1 (see Goldreich & Tremaine 1980; Ward 1988; Masset 2008; see also the appendix of Duffell & Chiang 2015 for a summary). These are at the location of the planet, that is, Ω = Ω_pl, and their strength decreases with m. The main resonances that drive the eccentricity are the ‘external’ first-order Lindblad resonances. These are located interior (for l = m + 1, $\text{Ω=}\left(\frac{m+1}{m-1}\right){\text{Ω}}_{\text{pl}}$) or exterior (for l = m − 1, $\text{Ω=}\left(\frac{m-1}{m+1}\right){\text{Ω}}_{\text{pl}}$) to the planet’s orbit, become closer and closer to the planet’s orbit for increasing m, and their sizes also initially increase with m, but ultimately become smaller for higher ms.

For small bodies (which do not open a gap), Artymowicz 1993 found that the e-damping due to the coorbital Lindblad resonances overcomes any excitation due to the external resonances by roughly a factor of three; instead, for massive enough planets, their gap is significantly deep at the corotation radius that the e-damping coorbital resonance’s effect is quenched. Sufficiently massive planets can even have their eccentricities even excited (Papaloizou et al. 2001; Kley & Dirksen 2006; Bitsch et al. 2013) due to the combined effect of the remaining resonances (e.g. Duffell & Chiang 2015). Thus, as an intermediate-mass planet carves a partial gap, we expect to transition smoothly from one regime (e-damping) to the other (e-pumping).

One might imagine that, as a gap is cleared at the location of the planet, the efficiency of e-damping resonances decreases in proportion to the decrease of gas surface density, thus explaining the decrease in e-damping efficiency. However, works on the shape and width of gaps (Crida et al. 2006; Kanagawa et al. 2018, as well as our own simulations) show that as soon as a gap opens, it already has a finite width, which remains relatively constant as the gap’s depth increases. Thus, as material is removed at the location of the e-damping resonances (i.e. the corotation radius) so is the case for some e-pumping resonances (the ‘external’ Lindblad resonances closest to the planet), so the picture is not immediately clear. However, we computed the effect of each specific first-order Lindblad resonance using disc profiles from our hydrodynamical simulation; we find that this effect is almost completely balanced by the overdensity at the inner and outer edges of the gap, while at the same time the locations closer to the planets are associated with higher indices m, which did not contribute much to the final ė. Thus, the resulting e-pumping efficiency remains relatively constant and decreases by less than 10%. Since the main term that governs ė/e is given by the e-damping terms corresponding to the co-orbital first-order Lindblad Resonances with l = m + 1, m − 1, and these are proportional to the gap depth in surface density, this explains why, down to a decrease of e-damping efficiency by roughly a factor of a few due to the carving of a partial gap at the location of the planet, we expect this decrease to be linear with Σ_min/Σ₀.

Appendix C Analysis of indirect forces

In our setup, the planets’ orbits are fixed and the forces from the disc to the planet are recorded. However, this situation is not purely consistent, since slight disc asymmetries can cause a non-vanishing acceleration onto the star, meaning that the reference frame with the star at the origin is non-inertial. These indirect forces might thus have an effect on the osculating elements, and therefore also on the eccentricity damping rates. In order to assess the effect of indirect forces, we track the acceleration onto the star due to its interaction with the disc. Since again the star is fixed onto the origin in the code’s reference frame, this results in an equal but opposite acceleration onto the planet. We see, however, that this acceleration is, on average, of the order of about two orders of magnitude smaller than the acceleration due to direct effects between the disc and the planet. We thus conclude that the asymmetries of the disc generated by the planet’s effect onto the disc itself are not large enough to cause large indirect forces’ contributions to the planet’s evolution. The picture would change in the case of more massive planets and more pronounced disc asymmetries.

Appendix D Fits to e-damping efficiency

In this appendix, we report and additional analysis of the data described in Section 4.3. Figure A.3(a) shows e-dependent (coloured dashed lines) and global (black dashed line) linear fits to the observed e-damping efficiency along the full gap-depth axis. While for small values of e/h this yields a good fit to the data, we see that this is not so for higher eccentricities, e/h ~ 1. The global fit (black dashed line) fails to accurately reproduce the data within acceptable errors for type-I migration formulas. For this reason, in Section 4.3 we resort to a double-linear fit that is eccentricity dependent. The left-part of the fit is shown in Figure A.3(b) for clarity. Here, the data used to make the fit have been truncated for observed gap depths of at least 80%. In this case, the combined fit alone (black dashed line) is able to fully recover the data within reasonable variations of less than 10%. The fit in Section 4.3 represents an extension of this approach, where an e-dependent fit to the right hand side of the plot (gap depths shallower than 80%) is continuously connected to the simple e-independent fit for deeper gaps.

References

All Tables

All Figures

Fig. 1 Comparison of orbital elements damping efficiencies 1/τmig (upper panel), 1/τe (middle panel), and 1/τa (lower panel) obtained from the output of our hydrodynamical simulations and from analytical formulas in the literature (CN2008 for Cresswell & Nelson 2008, P2011 for Paardekooper et al. 2010, and JM2017 for Jiménez & Masset 2017). The figures here are for the case of a mpl = 1 × 10−5M* planet with α = 3 × 10−4 and h = 0.05. τmig is obtained from the torque felt by the planet from the disc and typically shows good agreement between hydro and analytical predictions within the latter’s margin of uncertainty (20%). The analytical expression for τe comes from Cresswell & Nelson (2008) or from our fit presented in Sect. 4.3. τa is obtained from τmig and τe using Eq. (19).

In the text

Fig. 2 Predicted gap depth, Σmin/Σ0, from Kanagawa et al. (2018) versus value obtained from the simulations. Different colours represent different levels of viscosity as shown in the legend. The lower the viscosity, the bigger the difference in gap depth. However, down to the viscosities considered by Kanagawa et al. (2018), we see good agreement, as expected (the data for αt = 3.16 × 10−3 are shown with smaller symbols as they are only used for the analysis of the gap depth). Different shapes are used to represent different aspect ratios and planetary masses; squares, downward-pointing triangles, and upward-pointing triangles represent aspect ratios of 0.04, 0.05, and 0.06, respectively; empty, filled and crossed symbols represent mpl/M* = 10−5, 3 × 10−5, and 6 × 10−5, respectively.

In the text

Fig. 3 Surface density profile, Σt/Σ0, after ${\overline{t}}_{\max }=2950$ orbits and tmax = 3000 orbits of a mpl/M* = 3×10−5 planet in a disc with αt = 3.16 × 10−4 and h = 0.05. We note that after the last 50 orbits, the surface density has reached a steady state, with a change of less than 0.1%. The black horizontal dashed line represents the prediction of the gap depth from Kanagawa et al. (2018). The two dot-dashed horizontal lines are used here to measure the observed depth of the gap carved by the planet by taking the average of their values.

In the text

Fig. 4 Emergence of vortices in our simulations as a function of the gap depth. When no vortex was observed, we report both the observed gap depth at the end of our simulations, as well as the predicted one from Kanagawa et al. (2018), with a green circle (filled and unfilled, respectively). A red ⊗ symbolises those simulations where a vortex appeared, in which case only the predicted value from Kanagawa et al. (2018) is plotted. A dashed vertical line marks the approximate location of the transition between no vortex and vortex.

In the text

Fig. 5 Eccentricity damping efficiency versus gap-depth for all setups where no vortex emerged. In both panels, the e-damping efficiency on the vertical axis is normalised by the expected value from Cresswell & Nelson (2008). A 20% error around the expected value in the limit of no gap (to the right in the plots) is shown by two dashed horizontal grey lines. Values for the e-damping efficiency are shown for different eccentricities, e/h ∈ {0.2, 0.4, 0.6, 0.8, 1), by points of different colours joined together by opaque lines. Panel a shows on the horizontal axis the gap depth according to the prediction of Kanagawa et al. (2018), which results in an extremely noisy scatter plot. Panel b shows, on the horizontal axis, the gap depth observed from the actual simulations as plotted in Fig. 2, which shows a much cleaner dependence. In both panels, the overall trend clearly shows a decrease in e-damping efficiency (i.e. longer e-damping timescales) for deeper and deeper gaps, down to a factor of ~l/4 less efficient eccentricity damping at the transition from type-I to type-II regimes (gap depths of ≃0.25) as compared to the limit of no gap.

In the text

	Fig. 6 Double linear fit of e-damping efficiency to the data from panel b of Fig. 5 according to Eq. (28). The fit to the data is split into a linear fit over gap depths ∈ [0.3, 0.8], and a separate linear fit over gap depths ∈ [0.8,1], which continuously joins with the previous fit. The slopes for both fits are given in the legend in the top left corner.
In the text

Fig. 7 Similar to Fig. 2, but using Eq. (31) as an approximation for the gap depth, showing better agreement with the data. As in Fig. 2, different colours represent different levels of viscosity as shown in the legend, and different symbols represent different aspect ratios and planetary masses: squares, downward-pointing triangles, and upward-pointing triangles for h = 0.04, 0.05, and 0.06 respectively; empty, filled, and crossed symbols for mpl/M* = 10−5, 3 × 10−5, and 6 × 10−5, respectively.

In the text

	Fig. 8 Similar to Fig. 6, but using Eq. (31) as an approximation for the gap depth on the horizontal axis. The fit to the e-damping efficiency obtained from hydrodynamical simulations remains good, within the typical variations of analytical formulas for type-I migration.
In the text

	Fig. A.1 Surface density evolution in a setup with αt = 10−3 and h = 0.05, over 3000 orbits of a mpl/M* = 1 × 10−5 planet on a (fixed) circular orbit at 1 AU. We show the surface density contrast (∑t(r, θ) − 〈∑t(r)〉θ)/〈∑t(r)〉θ on the (θ, r) plane, where 〈∑t(r)〉θ is the azimuthally averaged surface density at the distance r.
In the text

	Fig. A.2 Evolution of the torque felt by the planet in the setup presented in Figure A.1 as a function of time. The vertical dashed line represents timestamps for which the surface density is shown in Figure A.1. The torque stabilises quickly in this case, as the system reaches a long-lived steady state.
In the text

Fig. A.3 Linear fit of the e-damping efficiency to the data from panel b of Figure 5, that is, using the gap depths observed from the simulations on the horizontal axis. In panel a we fit the full data, both for the separate values of e/h (coloured dashed lines) and for the full data set combining all eccentricities (black dashed line). In panel b we fit the data in the same way, but considering gap depths of up to 80% (indicated by the dashed red vertical line). In both cases we report the observed slopes of the lines in the top left legend.

In the text