A Recipe for Eccentricity and Inclination Damping for Partial-gap Opening Planets in 3D Disks

In a previous paper, we showed that, like the migration speed, the eccentricity damping efficiency is modulated linearly by the depth of the partial gap a planet carves in the disk surface density profile, resulting in less efficient e-damping compared to the prescription commonly used in population synthesis works. Here, we extend our analysis to 3D, refining our e-damping formula and studying how the inclination damping efficiency is also affected. We perform high-resolution 3D locally isothermal hydrodynamical simulations of planets with varying masses embedded in disks with varying aspect ratios and viscosities. We extract the gap profile and orbital damping timescales for fixed eccentricities and inclinations up to the disk scale height. The limit in gap depths below which vortices appear, in the low-viscosity case, happens roughly at the transition between classical type-I and type-II migration regimes. The orbital damping timescales can be described by two linear trends with a break around gap depths ∼80% and with slopes and intercepts depending on the eccentricity and inclination. These trends are understood on physical grounds and are reproduced by simple fitting formulas whose error is within the typical uncertainty of type-I torque formulas. Thus, our recipes for the gap depth and orbital damping efficiencies yield a simple description for planet–disk interactions to use in N-body codes in the case of partial-gap opening planets that is consistent with high-resolution 3D hydrosimulations. Finally, we show examples of how our novel orbital damping prescription can affect the outcome of population synthesis experiments.


INTRODUCTION
One of the goals of a planet formation model is to predict, for a given system or for the full exoplanet population in a statistical sense, what planetary system we expect to form inside a disk with given physical properties (such as surface density, temperature and thickness profiles, a given level of turbulent viscosity, etc.) orbiting a given star (Ida & Lin 2008;Mordasini et al. 2009;Alibert et al. 2013;Alessi et al. 2017;Izidoro et al. 2017;Ndugu et al. 2018;Bitsch et al. 2019;Guilera et al. 2019;Izidoro et al. 2021;Emsenhuber et al. 2021;Savvidou & Bitsch 2023).A fundamental ingredient in such a model is the description of planet-disk interactions.A planet embedded in a disk modifies the disk's structure and evolution, and, in turn, this interaction causes the planet's orbit to modify its size, shape, and orientation around the host star.
The presence of a gap in the disk surface density at the location of the planet's orbit is one of the most evident signatures of planet-disk interactions.The existence of gaps in protoplanetary disks is now observationally well established (e.g. the ALMAbased DSHARP survey, Andrews et al. 2018;Huang et al. 2018 ; see also Segura-Cox et al. 2020, andBae et al. 2023 for a review).Although gaps cannot be linked directly to gap-carving planets as a general rule (because other physical processes in disks can alone explain these features (Béthune et al. 2017;Suriano et al. 2019;Riols et al. 2020;Cui & Bai 2022), and some putative forming planetary systems would actually be dynamically unstable if to each gap corresponded a planet, e.g.Tzouvanou et al. 2023), the presence of a gap in the gas disk is shown in some cases to correspond to the presence of young forming protoplanets (Keppler et al. 2018;Wagner et al. 2018;Haffert et al. 2019;Pinte et al. 2019;Izquierdo et al. 2022).
If the mere presence of a dip in the surface density can be an indication of the presence of a planet, the depth of the gap that the planet carves (the ratio between the minimum of the surface density profile around the location of the planet and the surface density for the same unperturbed disk) is a crucial parameter to quantitatively describe the interaction between a planet and its surrounding protoplanetary disk.On the one hand, when a strong enough perturbation has occurred, the planet can create a pressure bump outside of its orbit which prevents pebbles from drifting inwards (Paardekooper & Mellema 2006;Morbidelli & Nesvorny 2012;Lambrechts et al. 2014;Ataiee et al. 2018;Bitsch et al. 2018;Weber et al. 2018); the planet mass at which this occurs is called the pebble isolation mass.In our Solar System, the fact that Jupiter's core would have stopped the inflow of pebbles from the outer disk (the so-called Jupiter barrier) could explain the observed isotopic dichotomy between non-carbonaceous and carbonaceous meteorites (Kruijer et al. 2017(Kruijer et al. , 2020)).The velocity deviation from a pure Keplerian rotation due to the acceleration of gas in response to the presence of a planet (so-called velocity kinks) can also be used to detect the presence of planets (e.g.Teague et al. 2018;Pinte et al. 2018Pinte et al. , 2019;;Izquierdo et al. 2022;Pinte et al. 2023).On the other hand, the carving of a gap and the establishment of a pressure bump affect in turn the planet's own physical and dynamical evolution.By stopping pebbles outside of its orbit it cannot efficiently accrete pebbles anymore so its solid budget is confined.Moreover, the local change in surface density (the gap) modulates the strength of planet-disk interactions that drive a growing planet's disk-driven migration.Classically, migration regimes are divided into a type-I regime, when the planet's mass is low enough (typically up to a few ten's of an Earth's mass for a typical protoplanetary disk) that it does not modify the disk's underlying surface density much, and a type-II regime, where the planet starts carving a significant enough gap (Crida et al. 2006;Kanagawa et al. 2018;Robert et al. 2018), influenced also by the accretion of gas onto the growing planet (e.g.Crida & Bitsch 2017;Bergez-Casalou et al. 2020).
One thus needs to understand the formation of gaps on quantitative grounds.We note that the gap depth is known to depend not only on the planet's mass but on disk properties as well, namely its aspect ratio and turbulent viscosity (Crida et al. 2006;Kanagawa et al. 2018).In particular, in the presence of a planet of a given mass, thinner and less viscous disks will respond with a deeper gap than thicker and more viscous disks.In fact, the processes that were thought to drive turbulent viscosities, such as the Magneto-Rotational Instability (MRI, Balbus & Hawley 1991), are quenched in large portions of the disk midplanes where planets form, and the remaining hydro-instabilities generate viscosities that are at least an order of magnitude lower than expected (Lyra 2014;Pfeil & Klahr 2021;Barranco et al. 2018, see Lesur et al. 2023 for a review).The analysis of observed disks also shows that results are best reproduced for similar low viscosities (see for example HL Tau and Oph163131 respectively by Pinte et al. 2016 andVillenave et al. 2022).Thus, even a low-mass planet that would traditionally be considered in the type-I regime may start opening a partial gap.For this reason, as convenient as the separation of migration regimes may be, we must be able to accurately describe modes of planetary-disk interactions that lie in between the two classical extremes.We note that planets that would fall in such transitional regimes have masses of a few to a few ten's of Earth's mass (so called Mini-Neptunes or Super-Earths, depending on whether they feature a thin gaseous atmosphere or not).Such planets do not exist in our own Solar System but appear to be the most common type of exoplanet in the galaxy (Mayor et al. 2011;Fressin et al. 2013;Petigura et al. 2013;Winn & Fabrycky 2015;Zhu et al. 2018;He et al. 2021;Lissauer et al. 2023), and represent the cores of giant planets in the core accretion model (Pollack et al. 1996).
Planet-disk interactions can be understood as a combination of changes in the planet's orbit's size, shape and inclination with respect to the disk mid-plane; these are usually referred to as a planet's migration (which is typically inward, thus associated to a damping of the planet's semi-major axis1 ), eccentricity damping and inclination damping, respectively.Concerning the migration efficiency, Kanagawa et al. (2018) has shown that, for circular and coplanar orbits, the transition between migration speeds in the classical type-I and type-II regime is modulated linearly by the depth of the gap carved by the planet.In a previous paper (Pichierri et al. 2023), we showed using high-resolution 2D locally isothermal hydrodynamical simulations that, for non-inclined planets, eccentricity damping efficiencies follow a similar trend, with a linear dependence on the gap depth (whose slope and intercept depend on the eccentricity).This fact is supported by theoretical grounds based on the gap profile opened by the planet and our understanding of the so-called eccentricity waves responsible for driving the eccentricity evolution of a planet embedded in a disk (Goldreich & Tremaine 1980;Tanaka & Ward 2004;Ward 1988;Masset 2008;Duffell & Chiang 2015).We found that e-damping efficiencies can be significantly lower than in the case of shallow gaps (Tanaka & Ward 2004;Cresswell & Nelson 2008); this finding bridges the gap between the classical regime of eccentricity damping that is typically associated to low-mass planets and the eccentricity pumping that is observed for high-mass planets in the type-II regime (Papaloizou et al. 2001;Kley & Dirksen 2006;Bitsch et al. 2013;Duffell & Chiang 2015).We also obtained a fitting formula to predict the gap depth of a planet of a given mass embedded in a disk with a given aspect ratio and viscosity that is formally similar to the one from Kanagawa et al. (2018), but gives a more accurate prediction for partial-gap opening planets in the low-viscosity regime, not probed by Kanagawa et al. (2018).In this paper, we extend our analysis to inclination damping efficiencies by using high-resolution 3D locally isothermal hydrodynamical simulations, and allowing the planet to lie on inclined orbits with respect to the disk mid-plane.At the same time, we refine our results on the eccentricity damping since there are known differences in e-damping efficiencies between 2D and 3D (Tanaka & Ward 2004).The fitting formula for the gap depth for partial-gap opening planets in low-viscosity disks is also investigated in the 3D case.Finally, we use the resulting gas surface density profiles obtained in our simulations to characterise the final mass that such planets would reach by accreting pebbles in their protoplanetary disk (the pebble isolation mass) and compare the outcome with known results from the literature (Bitsch et al. 2018).
The paper is organised as follows.Section 2 describes our disk model and the setups used in our hydro-dynamical simulations.Section 3 gives a comprehensive description of planet-disk interaction schemes in hydro-and N-body simulations: this includes how orbital damping timescales are extracted from hydro simulations and how they are implemented in N-body codes, including the transition from the type-I to the type-II migration regime.Section 4 describes the main results of our hydro simulations, yielding simple fitting formulas for the (partial) gap opened by a planet, and for the eccentricity and inclination damping timescales when planets open partial gaps in their surrounding disks, as a function of the orbital eccentricity and inclination.In Section 5 we compare our work to previous results on pebble isolation mass scaling laws, to yield a complete description of the dynamical interactions of partial-gap opening planets with their surrounding protoplanetary disks, and we discuss the implications of our orbital damping formulas in the context of population synthesis models.Finally, Section 6 summarises our results.
We then add a planet at a distance of 1 AU on a fixed orbit (we show in Appendix A that the results on orbital damping efficiencies obtained with a planet on a fixed orbit are comparable to the case of a planet let free to evolve in the disc).The mass of the planet increases smoothly from an initial value of 0 to its final masses (which varies across simulations, see below) over the course of 50 orbits.We consider planetary masses that are always below the thermal mass m th = h 3 M * , so that the disk-planet tidal perturbation does not drive local nonlinear shocks and can be treated linearly (Lin & Papaloizou 1986).We thus have surface density and temperature profiles fixed across all simulations, while α t , h and m pl /M * are left as free parameters.We consider α t values between 3.16 × 10 −5 and 10 −3 , since in the MRI-dead zone a residual turbulent viscosities can arise by purely hydrodynamical instabilities, with α t of order 10 −4 (e.g.Pfeil & Klahr 2021;Flock et al. 2020;Lesur et al. 2023), and observational constraints determine α t in disks to range between 10 −5 and 10 −2 (Pinte et al. 2016;Rafikov 2017;Dullemond et al. 2018;Flaherty et al. 2017Flaherty et al. , 2018;;Villenave et al. 2022).We take the disk's aspect ratio h ∈ {0.04, 0.05, 0.06} and planet masses in the Super-Earth/Mini-Neptune range, m pl /M * ∈ {1 × 10 −5 , 3 × 10 −5 , 6 × 10 −5 }.Compared to the 2D case, for a given disk structure and planetary mass we need to vary not only e/h but also i/h (where e is the planet's eccentricity and i is the inclination); moreover, 3D simulations are more costly than 2D ones.For this reason, in order to spare computational resources, we do not investigate a full grid of parameters {m pl /M * , α t , h} like we did in Pichierri et al. (2023), but we only consider setups which allow us to reach different gap depths at relatively equally spaced intervals to obtain a 3D version of the results obtained in Pichierri et al. (2023).
We give in Table 1 the list of all our disk setups.For each of the setups in the top entries, we varied the values of e/h and i/h independently in {0, 0.25, 0.5, 1}.We also run simulations for additional setups in the non-eccentric and non-inclined cases (bottom entries of Table 1), which we used to update our prediction for the gap depth.These total to 167 high-resolution 3D simulations.We do not consider higher values for the eccentricity and inclination for similar reasons as our 2D paper: the analytical formulas that we wish to compare our results to start breaking down at higher e/h and i/h; small single planets in such mass range have their eccentricities/inclinations damped by the disk so they are not expected alone to reach large e or i (Bitsch & Kley 2010, 2011;Cresswell & Nelson 2008); when multiple planets interact in a disk, e.g. by capturing in resonance via convergent migration, the expected capture eccentricities are of order h (Papaloizou & Szuszkiewicz 2005;Crida et al. 2008;Goldreich & Schlichting 2014;Deck & Batygin 2015;Pichierri et al. 2018), and the inclinations can be excited only by second order effects.Along the simulation, the planet's orbit (semi-major axis a, eccentricity e and inclination i) is kept fixed, as well as its mass m pl (except for the initial mass taper).
We note that in our simulations the reference frame is centered on the star and indirect forces should be considered.Similarly to our 2D study, we 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 disk; the disk feels indirect forces from the planet; the indirect forces of the disk onto itself is not included since we do not consider the disk's self-gravity.

Parameters
Values Like in our previous paper Pichierri et al. (2023), we run our numerical experiments using the fargOCA code (fargo with Colatitude Added; Lega et al. 2014) 2 , which is a 3D extension of the fargo code (Masset 2000), 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 used 512 grid cells with arithmetic spacing in radius for a disk extending from 0.5 to 2 AU, 2000 cells in azimuth for the full (0, 2π), and 70 cells in zenith for a disk with colatitude of 83 • .This gives square-ish cells at the location of the planet with δr ≃ δϕ ≃ δθ ≃ 0.003.This resolution is similar to our 2D runs, which we achieved by considering a slightly narrower radial domain in order to maintain a manageable total number of grid cells.Our disk is however still larger compared to the ones used in Jiménez & Masset (2017), while at the same time achieving a higher resolution.This ensures that, even for the smallest planetary masses, 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 performed resolution convergence tests as described in Appendix B. We used a smoothing length for the potential of the planet of r sm = s sm R H with s sm = 0.5.Finally, we used radially evanescent and vertically reflecting boundary conditions (de Val-Borro et al. 2006).

Orbital elements damping timescales
At regular time intervals, the fargOCA code outputs the (direct) gravitational force felt by the planet from the disk and the force felt by the star from the disk (indirect term).This indirect force emerges because of the asymmetry in the gas density resulting from the gas' response to the presence of the planet.Since our simulations are performed in a non-inertial astrocentric frame of reference, this force will result in an indirect (fictitious) response force felt by the planet, so we need to add it to the direct force felt by the planet from the gas.The resulting force F disk→pl describes the sum of direct and indirect planet-disk interactions, and thus the true force felt by the planet in an inertial reference frame.As in Pichierri et al. (2023), we use orbit-averaged forces, where the average is done over 20 points along the planet's orbit.
Following Burns (1976), we decompose the force F disk→pl into three components: where e R , e T and e N represent an orthonormal vector triad such that e R is in the direction of r pl , e T is inside the orbital plane and transverse to r pl , and e N is perpendicular to the orbital plane in the direction e R × e T , which is also the direction of the (orbital) angular momentum vector L = mr pl × v pl .Here, r pl is the position of the planet, v pl := ṙpl is its velocity, and m = (m pl M * )/(M * + m pl ) ≈ m pl is the reduced mass of the planet.We also introduce L := ∥L∥ the norm of the angular momentum vector, given by where µ = 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.The norm of the angular momentum is independent of the orbit's inclination i, which only dictates the orientation of L (see below).Finally, we introduce the (orbital) energy We need to determine how the different components of the perturbing force (1) translate into orbital elements damping.This damping is typically defined through damping timescales τ a , τ e and τ i via like in the 2D case we also define the migration timescale τ m via that is, the damping timescale of the (norm of the) angular momentum, which is different from the semi-major axis evolution timescale.Only forces inside the orbital plane (Re R + T e T ) can change the orbit's shape (semi-major axis and eccentricity, and thus the norm of the angular momentum), while forces perpendicular to it will change its orientation, i.e. the orbit's inclination.Inside the orbital plane, the evolution of a and e due to disk-planet perturbative forces is similar to the 2D case in Pichierri et al. 2023.The time derivatives of L = Le N and of its norm L are given by the torque: where r pl = ∥r pl ∥ and only r pl T e N contributes to the change in L since −r pl Ne T is perpendicular to L and only changes its direction (see below).We then consider the power which represents a change in orbital energy, P = Ė.Using these expressions, one calculates (see e.g.Pichierri et al. 2023 for an explicit derivation) where and C(e) ≈ e 2 for small eccentricities.For a circular orbit τ m = 2τ a .
To obtain τ i , we follow Burns (1976) equation ( 32) (see also Bitsch & Kley 2011) and write where θ pl is the planet's true longitude, and average this quantity over one orbit.Only N appears here because forces in the orbital plane cannot change the plane's orientation.This yields 3.2.3D planet-disk interactions in N-body codes N-body codes implement type-I migration using timescales τ m , τ e and τ i to define accelerations onto the planet given by (Papaloizou & Larwood 2000) where r pl and v pl are the planet's position and velocity, and k is the unit vector in the vertical direction.These accelerations directly damp canonical momenta (the three Delaunay action variables, e.g.Morbidelli 2002) associated to the orbital elements a, e and i as we show below.
The first equation describes a torque, i.e. a change in the angular momentum vector as L = mr pl × rpl = −m(r pl × v pl )/τ m ≡ −L/τ m .This implies that its norm L (the first Delaunay action) evolves according to (7).Since a m ∥ v pl , the resulting force lies on the orbital plane.
Equation ( 18) also represents a force that lies on the orbital plane, but it has zero torque since a e ∥ r pl , so it does not contribute to a change in L. By applying (18) over an orbit, the quantity E = (1 − e 2 ) −1/2 − 1 is damped exponentially over a timescale τ e /2 (Pichierri et al. 2023).E is the ratio between the Angular Momentum Deficit (AMD) Γ = m √ µa(1 − √ 1 − e 2 ) of the planet (Laskar 1997) and the norm of the angular momentum vector L, and it represents the second Delaunay action.Note that E ≃ e 2 /2 for small e's, so that Ė/E = −2/τ e translates into ė/e = −1/τ e for small e's.Therefore, at small e's, Equation ( 18) implements an exponential damping of the eccentricity over a timescale τ e as described by equation ( 5).The semi-major axis evolution described by (4) results from a combination of torque and e-damping with a timescale given by which is the equivalent of (13).Equation ( 19) implements a force such that Lx /L x = Ly /L y = −1/τ i and Lz /L z = 0, where L {x,y,z} are the component of L.
This quantity (up to a sign) represents the third Delaunay action.From the expressions of L{x,y,z} , one easily calculates that I follows the evolution İ/ For small I's, İ/I ≈ (−2/τ i ) so that I gets exponentially damped over a timescale τ i /2.Moreover, I ≃ i 2 /2 for small i's (i.e. for small I's), so İ/I ≈ −2/τ i translates to i/i ≈ −1/τ i for small i's.Therefore, at small i's, Equation ( 19) implements an exponential damping of the inclination over a timescale τ i as described by equation ( 6).
For inclined orbits, the acceleration a i in Equation ( 19) has components both inside and perpendicular to the orbital plane.For this reason, it also modifies the norm of the angular momentum, and therefore includes an additional unwanted (albeit small) change in the orbit's shape ( L 0).To overcome this nuisance, one can define a modified perturbing acceleration given by ãi where e N = (sin i sin Ω, − sin i cos Ω, cos i) ⊺ is again the versor orthogonal to the orbital plane and parallel to L. With this modified acceleration, Lx /L x = Ly /L y = −L z /(Lτ i ) (they gain a factor of cos i = L z /L), while L = 0 as desired.From this, one easily calculates that, under (21), the quantity I = 1 − cos i follows the evolution İ/I = (−1/τ i )(2 − I), which has the closed form solution For small I's, İ/I ≈ (−2/τ i ) so that also in this case I gets exponentially damped over a timescale τ i /2.Thus, like before, Equation ( 21) implements at small i's an exponential damping of the inclination over a timescale τ i as described by equation ( 6), but without introducing a spurious damping of L. This N-body implementation is thus better suited to be in line with the output of hydrodynamical simulations in the case of inclined orbits (at the expense of calculating the versor e N = L/L).
3.2.1.Type-I forces and transition to the type-II regime Analytical formulas for the damping timescales τ m , τ e and τ i for various migration regimes have been the subject of a number of works and we briefly summarise them here in the context of this work.
Much attention has been given to the expression of the torque Γ (i.e. of τ m by Eq. ( 11)) in the circular and non-inclined case for low-mass, non-gap opening planets (e.g., Tanaka et al. 2002;Cresswell & Nelson 2008;Paardekooper et al. 2011;Jiménez & Masset 2017).The total torque is typically split into a Lindblad component Γ L (usually negative) and a corotation component Γ C (usually positive, but prone to saturation).The total torque Γ tot,I = Γ L + Γ C provides a nominal type-I migration timescale τ m,I , which depends explicitly on the disk structure, particularly on the surface density profile (α Σ ), temperature profile (β T ), aspect ratio h and viscosity ν t at the location of the planet, as well as the planet's mass.The transition from non-gap opening (classic type-I) to partial gap opening planets to type-II planets has been investigated by Kanagawa et al. (2018), who found that the type-I migration timescales is modulated linearly by the depth of the gap carved by the planet in the disk: Here Σ min /Σ 0 measures the gap depth, where Σ min is the minumum of the (azimuthally averaged) surface density Σ(r) near the location of the planet r pl , while Σ 0 is the unperturbed disk surface density.Kanagawa et al. (2018) also provided a fitting formula for Σ min /Σ 0 : where is a dimensionless parameter.After the gap is considered to have been fully opened (Crida et al. 2006;Kanagawa et al. 2018), migration occurs in the type-II migration regime (Lin & Papaloizou 1986;Robert et al. 2018) and is beyond the scope of this work (we provide a measure of this limit in subsect.4.1).For eccentric and inclined orbits, the approach taken by many authors is to modulate the Lindblad and corotation torques by factors ∆ L and ∆ C that depend on e and i (Cossou et al. 2013;Pierens et al. 2013;Fendyke & Nelson 2014), to obtain a torque Γ tot,I = ∆ L Γ L + ∆ C Γ C for eccentric and inclined planets that can be used in population synthesis models (Izidoro et al. 2017(Izidoro et al. , 2021;;Emsenhuber et al. 2021).
The disk-driven eccentricity and inclination evolution has also been the subject of many works (e.g.Shu et al. 1983;Ward 1988;Artymowicz 1994;Ward & Hahn 1994;Tanaka & Ward 2004;Cresswell & Nelson 2008;Bitsch & Kley 2011;Pichierri et al. 2023).In particular, Tanaka & Ward (2004) gave the first expression for τ e and τ i , which was extended by Cresswell & Nelson (2008) to non-vanishing eccentricities and inclinations, by fitting the orbital evolution of a planet embedded in a 3D disk.Their fits yield Here τ wave is the typical type-I damping timescale (Tanaka & Ward 2004), were Ω K is the Keplerian orbital frequency and quantities with a subscript pl are evaluated at the position of the planet.These formulas are the most commonly used e-and i-damping prescription in the population synthesis literature (e.g.Izidoro et al. 2017, From Pichierri+23 Observed gap ○ Gap predicted by Pichierri+23 / Kanagawa+18 (no vortex in simulation) • Gap observed from simulation (no vortex in simulation) ⊗ Gap predicted from Pichierri+23 / Kanagawa+18 (vortex in simulation) Figure 1.Gap depth measured for the set of simulations of Table 1.When no vortex was observed, we mark both the observed gap depth and the predicted one with a green circle (filled and unfilled respectively).Predicted values are obtained with the prescription from Kanagawa et al. (2018) and from our 2D prescription from Pichierri et al. (2023), which, in the low viscosity case not probed by Kanagawa et al. (2018), gives a better prediction of the gap depth observed also in 3D simulations.In one simulation (for the setup marked with an asterisk in Table 1) we observe a vortex.This case is marked with a red circle with a cross.In this case we can not measure a gap depth (because the gap continuously changes); therefore, in the observed gap line, we simply use the value predicted in Pichierri et al. (2023).A dashed vertical line marks the approximate location where the transition between no vortex and vortex lies, which is similar to the one observed in 2D simulations (Pichierri et al. 2023).
-   26)).However, they are only valid for non-gap opening planets (Σ min /Σ 0 ≃ 1).In our preliminary 2D investigation (Pichierri et al. 2023), we extended the formula for τ e to partial gap opening, eccentric but non-inclined planets, and found that, like τ m , also τ e is modulated by a linear function of Σ min /Σ 0 , whose slope and intercept depend on e/h.In this work, we further extend this study to 3D simulations where planets are allowed to reside on inclined orbits too.

Gap opening and emergence of vortices
Like in our 2D experiments, we run all e = 0, i = 0 simulations up to at least 3000 planetary orbits (the integration time used in Kanagawa et al. 2018) and recorded the final gap depth Σ min /Σ 0 carved by the planet. 3We note that the gap will clear over a timescale 5 to 10 times (2/3)T pl (a pl /x s ), where T pl is the orbital period and x s is the half-width of the horseshoe region (Masset 2008), which is at most ∼ 2000 orbits in the cases we consider.The overdense regions located on each side of the gap will then spread radially to achieve a steady state, which is expected to happen over a longer (viscous) timescale.This timescale is too long to cover entirely in the context of hydro-dynamical simulations, so we adopt a practical approach where we prolong the integrations until the surface density does not vary significantly over time.Thus, in the lowest viscosity / thinnest disk runs, we extend the zero eccentricity and inclination run by an additional 1000 orbits.This ensures that, at the end of the simulations, the surface density always changes by less than 0.1% over 50 orbits, meaning they have reached a quasi-steady state.
In Pichierri et al. (2023), we found that an observed Σ min /Σ 0 ≃ 0.25 marked the transition from stable runs (Σ min /Σ 0 ≳ 0.25) and the appearance of vortices (Σ min /Σ 0 ≲ 0.25).We confirmed that a similar behaviour is observed in 3D simulations.Indeed, we run a setup with m pl /M * = 6 × 10 −5 , α t = 3.16 × 10 −5 and h = 0.05, which has a predicted gap depth of Σ min /Σ 0 ≃ 0.26, which resulted in an unstable disk (marked with an asterisk in Table 1 -we note that for such a setup, all runs with different eccentricities and inclinations became unstable); instead, when m pl /M * = 6 × 10 −5 , α t = 3.16 × 10 −4 and h = 0.04, which has a predicted gap depth of Σ min /Σ 0 ≃ 0.3 (which is also the observed gap depth), the disk remained stable.We note that, like in the 2D case, the vortex appears for the lowest viscosity and the most massive planet, when the gap depth approaches ∼ 0.25.Therefore the transition from type I migration to type II is also a transition from the no-vortex case to the vortex one but only because we are in a very low viscosity context.Figure 1 shows a diagram where we label the outcome of each simulation 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) or Pichierri et al. (2023); 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 Pichierri et al. (2023).Our 2D-fitting formula from Pichierri et al. ( 2023) is similar to Kanagawa et al. (2018)'s formula (23) (which was obtained for α t = 10 −3 ), but provides a better fit to the data for lower viscosities, using a modified K factor given by K2D = 3.93q 2.3 h −6.14 α −0.66 t .We show in Figure 2 (panel a) that this formula gives a good approximation for the gap depths obtained from 3D simulations as well.An improved fit that better matches the gap opening in 3D simulations is given by where K3D = 28q 2.3 h −5.4 α −0.72 t . (29) Figure 2 panel (b) shows how this formula fits the data from 3D simulations.
When deeper gaps are carved by the planet and vortices appear, we cannot draw any definitive conclusion on the value of Σ min /Σ 0 or of the orbital damping timescales.However, like in Pichierri et al. (2023), we note that planets that carve such deeps gaps are already considered in the type-II regime (Kanagawa et al. 2018;Crida et al. 2006), which is outside the scope of this work.We also stress that we use the gap depth in the circular and non-inclined case as representative of the gap depth carved by planets on eccentric and/or inclined orbits.This is justified in the limit of our analysis as Hosseinbor et al. (2007) showed that the gap carved by an eccentric planet is almost identical to the one carved by a planet on a circular orbit if e < (m pl /(3M * )) 1/3 , which is always the case for the setups considered here (see also Bitsch & Kley 2010).More recently, Sánchez-Salcedo et al. (2023) showed that the gap is fairly independent of the eccentricity if e ≲ h.Finally, Bitsch & Kley (2011) find little dependence of the gap depth for inclined orbits.

Eccentricity damping efficiency for partial-gap opening planets
We show in Figure 3 the observed eccentricity damping efficiency 1/τ e normalised by the expected efficiency 1/τ e,CN2008 from Cresswell & Nelson (2008) (Eq.( 25), which should be valid in the limit of no gap) as a function of the observed gap depth carved by the planet, and for different inclinations.Similarly to the 2D case, the data follow the following trends.For low eccentricities, 0 < e/h ≲ 0.5, one can fit 1/τ e as a function of the observed gap depth (in the limit of circular/non-inclined orbits) with a straight line over the full gap depth range of interest (Σ min /Σ 0 ≃ 0.3 to 1).In the limit of no gap (Σ min /Σ 0 ≃ 1), we recover the damping efficiency predicted by Cresswell & Nelson (2008) (itself based on the results of Tanaka & Ward 2004), as expected. 4For higher The different panels are for different orbital inclinations, with i/h ∈ {0, 0.25, 0.5, 1}.Two dashed horizontal gray lines indicate, around the expected value in the limit of no gap (to the right in the plots), an error of 20%, which is the typical uncertainty of analytical planet-disk interaction formulas (Paardekooper et al. 2011).For all inclination values, we observe a decrease in e-damping efficiency for deeper and deeper gaps well outside this margin of error, and down to a factor of ∼ 1/5 less efficient eccentricity damping at the transition from type-I to type-II regimes (gap depths of ≃ 0.3) as compared to the limit of no gap.The data are well modelled by a double linear fit that depends on the gap depth, the eccentricity and the inclination (see Eqs. ( 30) and ( 31)), shown with dashed lines whose color reflects the orbital eccentricity (the slopes for the piece-wise fits are given in the legend in the top left corner of each panel).
eccentricities, e/h ≃ 1, the data follow again a straight line for gap depths ≃ 0.3 up to ≃ 0.8, after which e-damping becomes significantly more efficient.This qualitative behaviour was already observed in Bitsch & Kley (2010); Fendyke & Nelson (2014) and in our 2D simulations Pichierri et al. (2023).The reason is that shallower gaps are also thinner in radial extent, and for sufficiently high e, the planet's excursions around r pl = a start interacting with the edge of the gap.Since the gap around a planet is carved where the Lindblad torques accumulate, around r pl ± 2/3H (Masset 2008), this happens when e/h ≃ 1.The specific linear dependence of the eccentricity damping efficiency as a function of the gap depth is also observed to depend on the orbit's inclination.Given the qualitative similarities with Pichierri et al. (2023), we obtain a fit to the data with a similar double-linear functional form which quantitatively reproduces the results of high-resolution 3D simulations:  The different panels are for different eccentricities, with e/h ∈ {0, 0.25, 0.5, 1}.Two dashed horizontal gray lines indicate, around the expected value in the limit of no gap (to the right in the plots), an error of 20%, which is the typical uncertainty of analytical planet-disk interaction formulas (Paardekooper et al. 2011).In all panels, there is a significant decrease in i-damping efficiency for deeper and deeper gaps, down to a factor of ∼ 1/5 less efficient damping at the transition from type-I to type-II regimes (gap depths of ≃ 0.3) as compared to the limit of no gap.The data are well modelled by a double linear fit that depends on the gap depth, the eccentricity and the inclination (see Eqs. ( 32) and ( 33)), shown with dashed lines of different colors depending on the orbital inclination.
where Figure 3 also shows a comparison between the fit (colored dashed lines) and the data (colored markers).This fit was obtained using the LinearModelFit function of the software package Mathematica to extract the linear models across different values of e and i as a function of the gap depth, and the NonlinearModelFit function to extract an explicit dependence on e and i. 5 The typical relative error given by the fit is of the order 10 -20% across all eccentricity and inclination values, which is of the order of the accuracy of torque formulas from the literature.We note that the slope of the fit in the vanishing eccentricity and inclination case is the same as our 2D study Pichierri et al. (2023).When Σ min /Σ 0 = 1 (no gap opened by the planet), the fit recovers known 3D eccentricity damping efficiencies exactly (Tanaka & Ward 2004;Cresswell & Nelson 2008).

Inclination damping efficiency for partial-gap opening planets
Repeating the same study for the inclination damping, we find that all the same arguments apply.Figure 4 shows the observed inclination damping efficiency 1/τ i normalised by the expected efficiency 1/τ i,CN2008 from Cresswell & Nelson (2008) (Eq.( 26)), as a function of the observed gap depth carved by the planet, and for different eccentricities.The overall trends are the same as in the case of the eccentricity damping, so we obtain a fit to the data with a similar functional form: where m2,i (e, i) = 0.8(e/h) + 1.12(e/h) 2 + 3.14(i/h − 0.25) − 0.42(e/h)(i/h − 0.25) + 2.9(i/h − 0.25) 2 . ( This double-linear fit is shown in Figure 4.The typical error given by the fit is less than 10% across all inclination and eccentricity values, and for Σ min /Σ 0 = 1 (no gap opened by the planet), we recover known 3D eccentricity damping efficiencies (Tanaka & Ward 2004;Cresswell & Nelson 2008).We note that, unlike in (31), there are coupling terms that are linear in e; this is in contrast with Cresswell & Nelson (2008)'s fit, which has no such terms.We checked that imposing that there be no linear coupling terms yields a worse fit to the data.Therefore, although we do not have a physical explanation for these terms, we use an agnostic approach and keep them in the fit for τ i to ensure the best possible match with the outcome of hydrodynamical simulations.

Pebble isolation mass
In our simulations, we set the planetary mass as a free parameter.However, in a real protoplanetary disk, this mass (at least in the limit of partial-gap opening planets where gas accretion is a negligeable effect) will be the result of accretion of solid, which is itself a process that depends on the gas structure and planet-disk interaction.In the pebble accretion scenario (Ormel & Klahr 2010;Lambrechts & Johansen 2012; see Johansen & Lambrechts 2017 for a review), the maximum mass that can be reached is the so-called pebble isolation mass (Morbidelli & Nesvorny 2012;Lambrechts et al. 2014;Ataiee et al. 2018;Bitsch et al. 2018;Weber et al. 2018).This is the mass at which the planet disturbs the disk surface density enough that it creates a pressure barrier outside of its orbit which prevents further pebbles to drift inwards and be accreted onto the planet.
Our high-resolution 3D simulations can also be used to validate previous works on the pebble isolation mass (Bitsch et al. 2018).From the output of our simulations, we define the gas density ρ(r) = Σ(r)/( √ 2πhr) and the pressure P(r) = c 2 s ρ(r), where c s = hrΩ K is the sound speed.We thus consider the pressure gradient ∂ log P ∂ log r and check when it changes sign in the vicinity of the planet's orbit.We then compare this outcome with the prediction from Bitsch et al. (2018) for the pebble isolation mass: We find that this prediction aligns well when our results, within an uncertainty of 20%, across all values of aspect ratios and viscosities considered here: when a planet is predicted to be below the pebble isolation mass, it does not generate a pressure barrier outside of its orbit in our simulation, and, conversely, when a planet is predicted to be above the pebble isolation mass, it generates a strong pressure bump.Note that Bitsch et al. (2018) also used 3D simulations run with the fargOCA code, albeit with a lower resolution.Thus, the prescriptions in Equations ( 30), ( 31), ( 32), ( 33) for the orbital damping timescales, together with Bitsch et al. (2018)'s formula (34) for the pebble isolation mass give a complete analytical description of the evolution of a super-Earth/Mini-Neptune planet embedded in a disk, from its limiting mass growth to its dynamical response to the presence of the disk, which is consistent with high-resolution 3D locally isothermal hydrodynamical simulations.2011)'s torque prescription; instead, panels on the left use Cresswell & Nelson (2008)'s e-damping prescription (eq.( 25)), while panels on the right use our modified e-damping prescription (eq.( 30)).The modified damping efficiency manifests itself in that, due to the overall less efficient e-damping, the captured resonant state is found at higher eccentricities; moreover, the libration inside the 3:2 resonance becomes overstable in the right panels, leading to an escape from the resonance and a more compact final state.

Application to N-body integrations
Analytical formulas for planet-disk interactions are widely used in the literature in the context of planet population synthesis models (Ida & Lin 2008;Mordasini et al. 2009;Alibert et al. 2013;Alessi et al. 2017;Izidoro et al. 2017;Ndugu et al. 2018;Bitsch et al. 2019;Izidoro et al. 2021;Emsenhuber et al. 2021).In particular, convergent migration (when two planets orbiting the same disk migrate in such a way that the sizes of their orbits approach each other) and eccentricity damping are generally associated to the assembly of mean motion resonant chains (Terquem & Papaloizou 2007;Cresswell & Nelson 2008;Morbidelli et al. 2008).However, the specific resonances that are built crucially depend on orbital damping efficiencies.This is important not only because the resonant structure attained at the end of the disk phase will be different on a quantitative level (i.e., which resonances will be observed in a given system), but also because the stability properties of these configurations are dependent on the resonance.More precisely, whether or not a given resonant chain assembled via disk-driven convergent migration will go unstable after the disappearance of the disk depends on how compact the chain is (Pichierri & Morbidelli 2020;Goldberg et al. 2022).Thus, although migration within a disk does in general lead to the assembly of resonant chains, which resonances are built has a strong impact in whether or not these resonance will even be observable after the the removal of the gas disk.
The processes of resonant capture and which resonances will be built under which conditions are fairly well understood (Batygin 2015; Deck & Batygin 2015;Pichierri et al. 2018;Batygin & Petit 2023).A mean motion resonance can be skipped if the resonance crossing time for the two planets is comparable to or shorter than the period of the planets' resonant interaction (the so called adiabatic limit, which is a condition on the torque), or if the dissipative torque is simply stronger than the resonant torque (which is a condition on the relative disk-driven e-damping onto the planets) (Batygin 2015;Batygin & Petit 2023).Moreover, even when the evolution is adiabatic and a resonance is successfully established, the presence of eccentricity damping breaks in general the adiabatic regime and the resonant equilibrium point may become an unstable fixed point, leading to so-called overstable librations and the escape from the resonant state (Goldreich & Schlichting 2014).For a fixed mass ratio between the planets, this is essentially a condition on the relative eccentricity damping efficiencies onto the two planets (Deck & Batygin 2015;Xu et al. 2018).Thus, even when the torques satisfy the adiabatic and stability conditions, different eccentricity damping efficiencies will lead to different final (resonant) states.
To elucidate this point, we show examples of N-body integrations of two planets inside a disk in Figures 5 and Figures 6 for the 3:2 and 4:3 resonance, respectively.These examples are not meant to depict realistic resonant capture scenarios, but rather to stress the differences that arise from the two eccentricity damping prescriptions.For this reason, we consider a constant disk surface density and aspect ratio, so that when the planets migrate inward their planet-disk interactions remain unchanged.We take α t = 3.16 × 10 −4 and h = 0.05.In Figure 5, we simulate two planets with masses m 1 /M * = 3 × 10 −5 and m 2 /M * = 5 × 10 −5 (which are below the pebble isolation mass for these disc parameters), starting with initially circular and coplanar orbits, and with initial period ratio slightly larger than 3:2.Since the surface density is constant at all radii, and planet 2 is more massive than planet 1, it will migrate inward faster than planet 1, so their period ratio will decrease and the two planets will approach the 3:2 commensurability under convergent migration.What happens after this is different in the panels on the right compared to the ones on the left.In the panels on the left, we implemented Paardekooper et al. (2011)'s torque prescription in conjunction with Cresswell & Nelson (2008)'s eccentricity damping formula, as is commonly done in the literature (e.g.Izidoro et al. 2017Izidoro et al. , 2021;;Emsenhuber et al. 2021).We observe that a successful capture has occurred (the evolution is in the adiabatic limit) and the resonant state achieved is stable (no overstable libration and jumping out of the resonance after capture).In the panels on the right, we used exactly the same initial conditions but we implemented our modified eccentricity damping formula (30).The planets have estimated gap depths of 0.87 and 0.67 respectively, so eccentricity damping is less efficient overall and the planets attain higher eccentricities (e-damping on planet 1 is very similar under both prescription, while planet 2 undergoes a less efficient e-damping with the modified formula (30), given its deeper gap).We observe that the final resonant state is overstable, that is, the amplitude of libration increases in time, and the system jumps out of the 3:2 resonance, to eventually end up in a more compact configuration.Figure 6 shows a similar result in the case of the 4:3 mean motion resonance, with planet masses m 1 /M * = 2×10 −5 and m 2 /M * = 5 × 10 −5 , where again α t = 3.16 × 10 −4 and h = 0.05.This case appears even more elusive at first, as the capture eccentricities appear to be very similar, but the stability properties of the resonant equilibrium are not.This is because in this case the eccentricity damping on the inner planet (which has a gap depth of 0.94) is enhanced, while that on the outer planet (which has a gap depth of 0.67) is less efficient using our modified formula (30) compared to the pure (Cresswell & Nelson 2008) prescription; thus, although the total eccentricity damping onto the planets is similar, the relative e-damping is different, the system undergoes overstable libration and exits the resonance.
In both cases, the final state will be more compact using the modified e-damping prescription (30) than using Cresswell & Nelson (2008)'s formula (25), leading to a system more susceptible to instabilities after the disk is removed (Pichierri & Morbidelli 2020;Goldberg et al. 2022).Even if the libration does not become overstable, the final eccentricities inside the same resonance may be higher for planets opening moderate gaps, which is also known to lead to less stable systems (Pichierri et al. 2018;Pichierri & Morbidelli 2020).Although simple by design, these experiments show that taking into account a more realistic modeling of orbital damping efficiencies for partial gap opening planets may have noticeable effects in the final product of population synthesis models.In particular, it may resolve the need to resort to more massive planets in order to trigger the instabilities needed to explain the orbital period distribution of known exoplanets (Izidoro et al. 2017(Izidoro et al. , 2021)).
Planetary inclination are affected by mean motion resonances only by second order effects (i.e., terms that are proportional to e × i 2 for first order resonances).In population synthesis models, they typically arise from seeding the planets with small initial inclinations, which may be subsequently excited by close encounters, collisions and scattering events.Because of the stochastic nature of this process, we do not systematically investigate the details of how our modified damping formulas (32), (33) might impact population synthesis calculations.In general, we expect that mutual inclinations will be enhanced, especially for planets opening deep gaps, because of the reduced i-damping efficiency in these cases.

CONCLUSIONS
In this paper we investigated planet-disk interactions for partial gap opening planets using high-resolution 3D hydro-dynamical simulations of locally isothermal disks with varying levels of turbulent viscosities and aspect ratios with an embedded planet of varying mass.The goals and methodology are similar to the ones used in a previous paper (Pichierri et al. 2023), which was limited to the 2D case: we reconsidered the problem of orbital damping timescales for planets that are classically in the type-I migration regime (a few to a few ten's of Earth's mass) but which would open partial gaps in disks of low viscosity and/or in thin disks, and are thus in between the type-I and type-II regimes.The expression of the torque felt by such planets in disks of arbitrary viscosities has been the subject of various works (e.g.Crida et al. 2006;Paardekooper et al. 2011;Jiménez & Masset 2017;Kanagawa et al. 2018) and these migration prescriptions have been used in population synthesis models to reproduce the observed characteristics of exoplanetary systems (e.g.Izidoro et al. 2017;Ndugu et al. 2018;Bitsch et al. 2019;Ogihara & Hori 2020;Izidoro et al. 2021;Emsenhuber et al. 2021).The transition between classical type-I and type-II torques for partial gap opening planets on circular and non-inclined orbits has been shown to depend linearly on the gap depth carved by the planet (Kanagawa et al. 2018), and we showed in Pichierri et al. (2023) that an equivalent linear trend is also observed in the eccentricity damping efficiency.In particular, e-damping efficiencies can be significantly lower (i.e., eccentricity damping timescales can be longer) than in the case of shallow gaps, which may have important consequences on the outcome of population synthesis models.This also fills the gap between the observed eccentricity damping that is typically associated to low-mass planets and the eccentricity pumping that is observed for very high mass planets (Papaloizou et al. 2001;Kley & Dirksen 2006;Bitsch et al. 2013).
Here, we extend the study to 3D disks and allow the planets to reside on orbits that are eccentric as well as inclined with respect to the disk mid-plane.We considered Super-Earth-type planets of varying fixed masses (m pl /M * = 1 to 6 × 10 −5 ) and varying orbital eccentricities and inclinations (e/h and i/h ranging from 0 to 1, where h = H/r is the disk's aspect ratio), embedded in disks of varying viscosities (3.16 × 10 −5 to α t = 10 −3 ) 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 period until a steady-state is achieved.
We analysed the surface density profile of the disk in response to the presence of the planet, in particular the depth of the gap opened by the planet and the threshold beyond which vortices appear, which gave similar results to the 2D case (Pichierri et al. 2023).Our fit for the gap depth agrees well with the one from Kanagawa et al. (2018) in the higher-viscosity disks, but it better reproduces the observed gap in the low-viscosity regime.We also considered the establishment of a pressure bump outside the orbit of the planet that would cause the inflow of pebbles to stop, thus halting the accretion of solid material; we found that our simulations agree well with previous results on the scaling of the pebble isolation mass with planetary and disk parameters (Bitsch et al. 2018).
We then considered the eccentricity and inclination damping efficiencies and their dependence on the gap depth.We found a similar qualitative behaviour as in our 2D study Pichierri et al. (2023).The orbital damping efficiencies (rescaled by the expected efficiencies in the no-gap case, Tanaka & Ward 2004;Cresswell & Nelson 2008) are well described as linear functions of the gap depth with slopes and intercepts that depend in general on the eccentricity and inclination; a break is observed around gap depths of 80% after which, for shallower gaps, the damping efficiency's slope with respect to the gap depth increases (see Fig. 's 3 and  4).These features can be understood on theoretical grounds (Pichierri et al. 2023).We therefore used an equivalent functional form as our 2D fit from Pichierri et al. (2023) and obtained an explicit but simple formula that depends on the gap depth, the orbital eccentricity and inclination (Eq.'s (30), ( 31) and ( 32), (33)), and which approximates the outcome of 3D high-resolution hydrodynamical simulations within the errors of torque formulas commonly used in population synthesis works.This gives a simple but complete description of planet-disk interactions for partial gap opening planets (from the traditional type-I regime down to gaps close to the traditional type-II regime) to be used in N-body simulations.
Finally, we tested the consequences of our novel formulas in the context of planet population synthesis simulations, in particular in the formation of mean motion resonant chains.Resonances are naturally associated with convergent migration (when two planets orbit the same star embedded in the same protoplanetary disk and the sizes of their orbits change in such a way that the orbits get closer to each other, Terquem & Papaloizou 2007;Cresswell & Nelson 2008;Morbidelli et al. 2008).This process is well understood on theoretical grounds (e.g., Batygin & Morbidelli 2013;Batygin 2015), and in particular it is known that which mean motion resonances are skipped and which ones are successfully established depends on the orbital damping timescales (Batygin 2015;Deck & Batygin 2015;Xu et al. 2018;Batygin & Petit 2023), and so do the final eccentricities after successfully capturing in a given mean motion resonance (Papaloizou & Szuszkiewicz 2005;Crida et al. 2008;Goldreich & Schlichting 2014;Deck & Batygin 2015;Pichierri et al. 2018).We show simple examples in which our modified formulas for orbital damping timescales for partial-gap opening planets yield dynamically different results than the prescriptions used so far in the literature, and we stress in what way this would impact the orbital states obtained at the end of population synthesis models.In particular, the establishment of more dynamically excited and compact states may resolve the necessity for more massive planets in order to trigger the instabilities that can explain the orbital period distribution of known exoplanets (Izidoro et al. 2017(Izidoro et al. , 2021)).

ACKNOWLEDGMENTS
The authors are grateful to the anonymous referee for comments which improved the clarity and content of the manuscript.G. P. and B. B. thank the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support.G. P. also thanks the Barr Foundation for their financial support, and K. Batygin for helpful comments that improved the manuscript.E. L. whishes to thank Alain Miniussi for the maintainance and re-factorization of the code fargOCA.We acknolewdge HPC resources from GENCI DARI n.A0140407233.

APPENDIX A. FIXED VS. FREE PLANETS
In this section we compare the eccentricity and inclination damping timescales that one would infer by i) keeping the planet on a fixed orbit, or ii) letting the planet respond to the disk and fitting the evolution of the osculating orbital elements.The first method is the one we used in this paper, and which is also typically used to obtain the strength of the torque (Paardekooper et al. 2011;Jiménez & Masset 2017); the advantage of this method is that one can wait arbitrarily long until a steady state is reached without having to worry about the planet's orbital state changing over time.The second method is the one used e.g. in Cresswell & Nelson (2008) and links more directly to practical applications for population synthesis works and N-body integrations.
To check if the two methods agree, we run additional hydro-dynamical simulations where we track the time evolution of the orbital elements after releasing a planet on an initially eccentric and/or inclined orbit.In particular, we use as initial conditions (for both the gas and the planet) the end state of our simulations where the planet had been kept on a fixed orbit, after a steadystate has been reached, and we release the planet.We then set up N-body integrations mimicking planet-disk interactions with the same initial conditions as the hydro-simulations and check whether the evolutions of the eccentricity/inclination over time match with that of the hydro-simulations, that is, whether they are damped on a timescale comparable to the expected one.In these N-body simulations, we use both the classical Cresswell & Nelson (2008), and our modified prescription from Eq.'s (30), ( 31) and ( 32), (33) to mimic disk-driven e-and i-damping.
Since the typical type-I damping timescale τ wave is inversely proportional to m pl (Eq.( 27)), in order to spare computational resources we consider the cases with the highest planetary mass m pl /M * = 6 × 10 −5 .We then first consider disk parameters such that our fitting formulas do not constitute a significant deviation from Cresswell & Nelson (2008) (e.g.h = 0.05, α t = 10 −3 ) in order to make a fair comparison.Secondly, we consider one case where we instead expect there to be a noticeable difference in damping timescales (e.g.h = 0.04, α t = 3.16 × 10 −4 , which gives an estimated gap depth of ≃ 0.3, and thus a factor ∼ 4 difference in damping efficiencies).We consider as initial conditions e/h and/or i/h equal to 1, and we run the hydrodynamical simulations for 20 orbits, which is enough to track the damping in the eccentricities/inclinations.The results of the m pl /M * = 6 × 10 −5 , h = 0.05, α t = 10 −3 setup are presented in Figure 7, which shows the evolution of the semi-major axis and eccentricity and/or inclination for a planet with different combinations of initial conditions (eccentric but coplanar with the disk, inclined but circular, eccentric and inclined).In this case, as expected, the outcomes of N-body integrations with the standard Cresswell & Nelson (2008) prescription (labeled "CN2008, orig.") and with the modified prescription from this paper (labeled "CN2008,  .Evolution over time (in units of the orbital period T pl at a = a(0) = 1) of the orbital elements for a planet subject to disc-driven migration and e-and i-damping.The planet and disk parameters are reported above each panel.In all panels we plot the outcome of a hydrosimulation (labeled "fargOCA") and of two N-body integrations (one with the classical Cresswell & Nelson (2008) damping prescription, labeled "CN2008, orig.", and one with our modified prescription for τ e and τ i , labeled "CN2008, mod."), as shown in the legend at the bottom.Each column represents a different initial condition for the eccentricity and inclination of the planet, where either e, or i, or both, are initialised at h.We plot on the top row the evolution of the semi-major axis (notice that the curves for the two N-body integrations largely overlap) and on the bottom rows that of the eccentricity and/or inclination.Notice the log scale on the vertical axis, which we use in order to observe the slope of the curves, which is a measure of τ e and τ i ; the small oscillations in orbital elements (noticeable especially in the eccentricity evolution) have a frequency equal to the orbital frequency and have a relatively constant amplitude.
mod.") are rather similar to each other6 .They also match well against the outcome of the hydro-dynamical simulations with the released planet.When Cresswell & Nelson (2008)'s and our prescription differ more significantly (e.g., in the left and middle columns), the latter gives a better match to the outcome of our hydro-dynamical simulations.In a similar fashon, the results of the m pl /M * = 6 × 10 −5 , h = 0.04, α t = 3.16 × 10 −4 setup are shown in Figure 8.Here, the difference in evolution for the N-body integrations with the different damping prescriptions is more apparent as expected.We see that our modified prescription which takes into account the partial gap opened by the planet yields a good match to the outcome of these hydro-dynamical simulations.This analysis also shows that, within the parameter space considered in this work, whether a planet is kept on a fixed orbit or whether it is allowed to move in response to the gas would result in comparable orbital damping timescales.

B. CONVERGENCE TESTS
We run a few resolution tests to check the convergence of the results of our hydro-dynamical simulations.In all cases, the azimuthal extent is the full [0, 2π) interval and the colatitude is 83 • .
All three gave similar results.At low α ∼ 10 −4 , the lower resolution run appears slightly different from the higher resolution ones in that the disk goes mildly unstable, while Res.N and Res.H were still extremely similar.We thus use Res.N as our nominal resolution in order to spare computational resources while maintaining a good accuracy in our results.Figure 9 shows two examples of our resolution tests.

Figure 2 .
Figure 2. Panel (a): Gap depth prediction (either from Kanagawa et al. 2018 or from Pichierri et al. 2023, indicated by small markers linked together by dot-dashed lines or larger markers linked by continuous lines, respectively) versus the observed 3D gap depth.Panel (b): New 3D gap depth prediction (Eqs.(28), (29)) versus observed 3D gap depth in the limit of circular and non-inclined orbits.In all panels different colours represent different levels of viscosity according to the legend, and symbols of 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.

Figure 3 .
Figure3.Observed eccentricity damping efficiency versus observed gap-depth for all setups where no vortex emerged.In all panels, the e-damping efficiencies on the vertical axis are normalised by the expected value fromCresswell & Nelson (2008) and are shown for different eccentricities e/h ∈ {0.25, 0.5, 1} by points of different colors joined together by opaque lines.The different panels are for different orbital inclinations, with i/h ∈ {0, 0.25, 0.5, 1}.Two dashed horizontal gray lines indicate, around the expected value in the limit of no gap (to the right in the plots), an error of 20%, which is the typical uncertainty of analytical planet-disk interaction formulas(Paardekooper et al. 2011).For all inclination values, we observe a decrease in e-damping efficiency for deeper and deeper gaps well outside this margin of error, and down to a factor of ∼ 1/5 less efficient eccentricity damping at the transition from type-I to type-II regimes (gap depths of ≃ 0.3) as compared to the limit of no gap.The data are well modelled by a double linear fit that depends on the gap depth, the eccentricity and the inclination (see Eqs. (30) and (31)), shown with dashed lines whose color reflects the orbital eccentricity (the slopes for the piece-wise fits are given in the legend in the top left corner of each panel).

Figure 4 .
Figure 4. Similar to Figure 3, but showing the observed inclination damping efficiency (normalised by the expected value from Cresswell & Nelson 2008) versus the observed gap depth.Values for the i-damping efficiency are shown for different inclinations i/h ∈ {0.25, 0.5, 1} by points of different colors joined together by opaque lines.The different panels are for different eccentricities, with e/h ∈ {0, 0.25, 0.5, 1}.Two dashed horizontal gray lines indicate, around the expected value in the limit of no gap (to the right in the plots), an error of 20%, which is the typical uncertainty of analytical planet-disk interaction formulas(Paardekooper et al. 2011).In all panels, there is a significant decrease in i-damping efficiency for deeper and deeper gaps, down to a factor of ∼ 1/5 less efficient damping at the transition from type-I to type-II regimes (gap depths of ≃ 0.3) as compared to the limit of no gap.The data are well modelled by a double linear fit that depends on the gap depth, the eccentricity and the inclination (see Eqs. (32) and (33)), shown with dashed lines of different colors depending on the orbital inclination.

Figure 5 .
Figure 5. Evolution under convergent migration of two planets, with masses m 1 = 3 × 10 −5 M * and m 2 = 5 × 10 −5 M * respectively, near the 3:2 commensurability.The disk has α t = 3.16 × 10 −4 and h = 0.05.Panels on the left and on the right have the same initial conditions, and show, from top to bottom, the evolution of the semi-major axes, eccentricity, period ratio and resonant angles.All panels make use of Paardekooper et al. (2011)'s torque prescription; instead, panels on the left use Cresswell & Nelson (2008)'s e-damping prescription (eq.(25)), while panels on the right use our modified e-damping prescription (eq.(30)).The modified damping efficiency manifests itself in that, due to the overall less efficient e-damping, the captured resonant state is found at higher eccentricities; moreover, the libration inside the 3:2 resonance becomes overstable in the right panels, leading to an escape from the resonance and a more compact final state.
pl /M*=6×10 -5 , α=10 -3 , h=0.05 e pl (0)=0.05,i pl (0 Figure7.Evolution over time (in units of the orbital period T pl at a = a(0) = 1) of the orbital elements for a planet subject to disc-driven migration and e-and i-damping.The planet and disk parameters are reported above each panel.In all panels we plot the outcome of a hydrosimulation (labeled "fargOCA") and of two N-body integrations (one with the classicalCresswell & Nelson (2008) damping prescription, labeled "CN2008, orig.", and one with our modified prescription for τ e and τ i , labeled "CN2008, mod."), as shown in the legend at the bottom.Each column represents a different initial condition for the eccentricity and inclination of the planet, where either e, or i, or both, are initialised at h.We plot on the top row the evolution of the semi-major axis (notice that the curves for the two N-body integrations largely overlap) and on the bottom rows that of the eccentricity and/or inclination.Notice the log scale on the vertical axis, which we use in order to observe the slope of the curves, which is a measure of τ e and τ i ; the small oscillations in orbital elements (noticeable especially in the eccentricity evolution) have a frequency equal to the orbital frequency and have a relatively constant amplitude.

Figure 8 .ΓFigure 9 .
Figure8.Similar to Fig.7, but for planet and disk parameters such that the gap carved by the planet is deep enough that we would expect noticeable e-and i-damping efficiencies between the classicalCresswell & Nelson (2008) prescription and our modified damping formulas.Indeed, we observe less efficient damping in all cases, with our prescription following very closely the eccentricity and inclination evolution of the hydro-dynamical simulations.

Table 1 .
Set of parameters of our simulations.In the case marked with an asterisk, the disk was unstable.