A recipe for orbital eccentricity damping in the type-I regime for low viscosity 2D-discs

It is known that gap opening depends on the disc's viscosity; however, eccentricity damping formulas have only been derived at high viscosities, ignoring partial gap opening. We aim at obtaining a simple formula to model $e$-damping of the type-I regime in low viscosity discs, where even small planets may start opening partial. We perform high resolution 2D locally isothermal hydrodynamical simulations of planets with varying masses on fixed orbits in discs with varying aspect ratios and viscosities. We determine the torque and power felt by the planet to derive migration and eccentricity damping timescales. We first find a lower limit to the gap depths below which vortices appear; this happens roughly at the transition between type-I and type-II regimes. For the simulations that remain stable, we obtain a fit to the observed gap depth in the limit of vanishing eccentricities that is similar to the one currently used in the literature but is accurate down to $\alpha=3.16\times 10^{-5}$. We record the $e$-damping efficiency as a function of the observed gap depth and $e$: when the planet has opened a deep enough gap, a linear trend is observed independently of $e$; at shallower gaps this linear trend is preserved at low $e$, while it deviates to more efficient damping when $e$ is comparable to the disc's scale height. Both trends can be understood on theoretical grounds and are reproduced by a simple fitting formula. Our combined fits yield a simple recipe to implement type-I $e$-damping in $N$-body for partial gap opening planets that is consistent with high-resolution 2D hydro-simulations. The typical error of the fit is of the order of a few percent, and lower than the error of type-I torque formulas widely used in the literature. This will allow a more self-consistent treatment of planet-disc interactions of the type-I regime for population synthesis models at low viscosities.


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, andMAPS 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 planetdisc 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 A148, page 1 of 14 Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0),which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.This article is published in open access under the Subscribe to Open model.Open access funding provided by Max Planck Society.viscosity, etc; e.g.Ida & Lin 2008;Mordasini et al. 2009;Alibert et al. 2013;Alessi et al. 2017;Izidoro et al. 2017Izidoro et al. , 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(Izidoro et al. , 2021;;Pichierri & Morbidelli 2020;Goldberg et al. 2022), although a few remained stable in resonance and are still observed today 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. 2010Paardekooper et al. , 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. 2017Izidoro et al. , 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).Hydrodynamical 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 lowviscosity 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(Izidoro et al. , 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, A148, page 2 of 14 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 etc.).Moreover, 2D hydrodynamical simulations can be made to reproduce the main features of 3D simulations (Kley et al. 2009(Kley et al. , 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.

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 Σ(r) = Σ 0 (r/r 0 ) −α Σ and a constant aspect ratio of H/r = h = h 0 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 = α t H 2 Ω 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πα t h 2 r 2 Ω 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 (r) = T 0 (r/r 0 ) −β 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 edamping 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 −5 , 1. × 10 −4 , 3.16 × 10 −4 , 1. × 10 −3 }.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 −5 , 3 × 10 −5 , 6 × 10 −5 }, 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 3 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, 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) 1 .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 timeexplicit-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 0 = 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).

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 −3 ), 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 = 0.1, and derive a condition for gap opening given by 3 4 where R H is the Hill radius of the planet, q = m pl /M * , and R = r 2 pl Ω K,pl /ν t is the Reynolds number.Kanagawa et al. (2018) gave a prediction for the value of Σ min /Σ 0 , 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 where is a dimensionless parameter.Gyeol Yun et al. ( 2019) slightly improved this result, replacing Eq. ( 2) with Σ min /Σ 0 ≃ 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 −4 ); 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 −1 , 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 pl = 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, Σ t max (r)/Σ 0 (r), 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 Σ t max /Σ 0 (r) , 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 /Σ 0 .Our results are presented in Sect.4.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, and the power, 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, Γ := dL pl dt = Lpl , that is, the rate of change of the (orbital) angular momentum, which is given by where m = (m pl M * )/(M * + m pl ) ≈ m pl is the reduced mass of the planet, µ = G(M * + m pl ) ≈ GM * 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 := dE pl dt = Ėpl , that is, the rate of change of the (orbital) energy, which is given by A148, page 4 of 14 By the fundamental equation Γ = Lpl , the angular momentum of the planet evolves, and we can introduce a migration timescale τ m defined as Lpl 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 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: with C(e) ≈ e 2 for small eccentricities.Thus, we have We also see that for a circular orbit τ m = 2τ a .

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): where r pl and u 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, 2 an orbital damping of the eccentricity as described by Eq. ( 10).The semimajor axis evolution described by Eq. ( 9) thus results from a combination of torque and e-damping with a timescale given by that is, the equivalent of Eq. ( 16).

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 −3 .Defining the typical type-I damping timescale τ wave (Tanaka & Ward 2004), their best fit yielded (in the planar case) where is a reduction factor due to the planet's eccentricity and ê = e/h, while τ e,CN2008 = τ wave 0.780 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 2 By applying Eq. ( 18) over an orbit, the quantity that is damped exponentially (over a timescale τ e /2) is E = (1 − e 2 ) −1/2 − 1 ≃ e 2 /2 for small es.The quantity E is the ratio between the angular momentum deficit (AMD) of the planet (Laskar 1997) and the norm of the angular momentum vector.A148, page 5 of 14 A&A 670, A148 (2023) behaviour of the torque observed in their 2D numerical experiments with an error of up to 20%.Jiménez & Masset (2017) reanalysed 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(Izidoro et al. , 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), and a reduction of the corotation torque given by where e f = 0.5h + 0.01 is defined in Fendyke & Nelson (2014).The total torque is thus ).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.

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 ).The figures here are for the case of a m pl = 1 × 10 −5 M * 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).
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 /Σ 0 : 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 A148, page 6 of 14 G. Pichierri et al.: A recipe for orbital eccentricity damping in the type-I regime for low-viscosity 2D discs attempt to obtain an equivalent adaptation of Eq. ( 27) for the eccentricity damping for partial-gap opening planets.

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 −4 , m pl /M * = 6 × 10 −5 ), (h = 0.04, α t = 3.16 × 10 −5 , m pl /M * = 6 × 10 −5 ), (h = 0.04, α t = 10 −4 , m pl /M * = 6 × 10 −5 ), (h = 0.05, α t = 3.16 × 10 −5 , m pl /M * = 3 × 10 −5 ), and (h = 0.05, α t = 3.16 × 10 −5 , m pl /M * = 6 × 10 −5 ), 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.give this in the next section in terms of the depth of the gap carved by the planet.The case of (h = 0.05, α t = 10 −3 , m pl /M * = 1 × 10 −5 ), where a vortex does not appear, is shown as an example in Figs.A.1 and A.2.

Partial gap opening at low viscosity and/or thin discs
We calculated the gap depth, Σ min /Σ 0 , 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 −3 , while we go as low as α t = 3.16 × 10 −5 .In order to make a fairer comparison, we run additional hydrodynamical simulations with the same m pl and h ranges from Table 1 for α t = 3.16 × 10 −3 , 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 −3 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 −5 , α t = 3.16 × 10 −4 , and h = 0.05; the observed gap depth after 3000 orbits is ≃0.9, while the  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 m pl /M * = 10 −5 , 3 × 10 −5 , and 6 × 10 −5 , respectively.Σ t /Σ 0 Fig. 3. Surface density profile, Σ t /Σ 0 , after tmax = 2950 orbits and t max = 3000 orbits of a m pl /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 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.

Kanagawa+ gap depth prediction
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 A148, page 7 of 14 A&A 670, A148 ( 2023) Observed gap (when no vortex) ○ Gap predicted from Kanagawa+ (no vortex in simulation) • Gap observed from simulation (no vortex in simulation) ⊗ Gap predicted from Kanagawa+ (vortex in simulation) 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.
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 ≲ 0.25.

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. 2017Izidoro et al. , 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 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 ≃ 0.3 to 1).In the limit of no gap (Σ min /Σ 0 ≃ 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: 1 where m = max {m 1 , m 1 + 3(e/h − 0.3)}.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 /Σ 0 = 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.

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 Obs.gap depth 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 ∼1/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.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.
about 10 −3 , the prediction of Kanagawa et al. (2018) also shows relatively good agreement with our numerical experiments (with a perfect fit at α t = 10 −3 , 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 and we fit a function of the form K(q, α t , h) = Cq p q α p α t t h p h using our observed values for the right hand side of Eq. ( 29).We used   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 upwardpointing triangles for h = 0.04, 0.05, and 0.06 respectively; empty, filled, and crossed symbols for m pl /M * = 10 −5 , 3 × 10 −5 , and 6 × 10 −5 , respectively.
Mathematica's NonlinearModelFit function and obtained the following fit: K = 3.93q 2.3 h −6.14 α −0.66 t . (30) We show, in Fig. 7, an analogue of Fig. 2, but using K from Eq. ( 30) instead of K from Eq. (3), showing much better agreement.This suggests a modification to the gap depth formula, namely 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.
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-tomoderate 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 −3 .Also, one should always keep in mind that these studies have been carried out using different hydrodynamical codes.

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 −1 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 hydrodynamical 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.

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(Izidoro et al. , 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 −3 ) 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 −4 (e.g.Pfeil & Klahr 2021;Flock et al. 2020), and observational constraints determine α t in discs to range between 10 −4 to a few 10 −2 (Rafikov 2017;Dullemond et al. 2018;Flaherty et al. 2017Flaherty et al. , 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 −3 (Hühn et al. 2021).
A148, page 10 of 14 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 −5 ) 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(Izidoro et al. , 2021)), allowing instabilities to also develop for planets of lower masses than the ones needed in these works to efficiently break the resonant chains.

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. 2017Izidoro et al. , 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 −3 ) 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 −4 ).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 highresolution, 2D, locally isothermal hydrodynamical simulations of super-Earth-type planets of varying masses (m pl /M * = 1 to 6 × 10 −5 ) 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 −3 down to 3.16 × 10 −5 ) 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 Σ t max /Σ 0 , where Σ t max is the radial surface density profile at the end of the simulation and Σ 0 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 −3 ), 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 gapdepth fit of Eq. ( 31), this allows us to provide a simple A148, page 11 of 14 A&A 670, A148 (2023) prescription for eccentricity-damping in the type-I regime for partial-gap opening planets that is consistent with highresolution 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.

)Figure 1
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 fromCresswell & Nelson (2008).We compile our results for all of our numerical sample in Sect.4.3.

Fig. 1 .
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 m pl = 1 × 10 −5 M * 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 fromCresswell & Nelson (2008) or from our fit presented in Sect.4.3.τ a is obtained from τ mig and τ e using Eq.(19).

Fig. 2 .
Fig.2.Predicted gap depth, Σ min /Σ 0 , fromKanagawa 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 byKanagawa 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 m pl /M * = 10 −5 , 3 × 10 −5 , and 6 × 10 −5 , respectively.
et al.: A recipe for orbital eccentricity damping in the type-I regime for low-viscosity 2D discs

Fig. 6 .
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 .
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 upwardpointing triangles for h = 0.04, 0.05, and 0.06 respectively; empty, filled, and crossed symbols for m pl /M * = 10 −5 , 3 × 10 −5 , and 6 × 10 −5 , respectively.

Table 1 .
Set of free parameters of our simulation.