Trapping (sub-)Neptunes similar to TOI-216b at the inner disk rim Implications for the disk viscosity and the Neptunian desert

Context. The occurrence rate of observed sub-Neptunes has a break at 0.1 au, which is often attributed to a migration trap at the inner rim of protoplanetary disks where a positive co-rotation torque prevents inward migration. Aims. We argue that conditions in inner disk regions are such that sub-Neptunes are likely to open gaps, lose the support of the co-rotation torque as their co-rotation regions become depleted, and the trapping e ﬃ ciency then becomes uncertain. We study what it takes to trap such gap-opening planets at the inner disk rim. Methods. We performed 2D locally isothermal and non-isothermal hydrodynamic simulations of planet migration. A viscosity tran- sition was introduced in the disk to (i) create a density drop and (ii) mimic the viscosity increase as the planet migrated from a dead zone towards a region with active magneto-rotational instability (MRI). We chose TOI-216b as a Neptune-like upper-limit test case, but we also explored di ﬀ erent planetary masses, both on ﬁxed and evolving orbits. Results. For planet-to-star mass ratios q (cid:39) (4–8) × 10 − 5 , the density drop at the disk rim becomes reshaped due to a gap opening and is often replaced with a small density bump centred on the planet’s co-rotation. Trapping is possible only if the bump retains enough gas mass and if the co-rotation region becomes azimuthally asymmetric, with an island of librating streamlines that accumulate a gas overdensity ahead of the planet. The overdensity exerts a positive torque that can counteract the negative torque of spiral arms. Under suitable conditions, the overdensity turns into a Rossby vortex. In our model, e ﬃ cient trapping depends on the α viscosity and its contrast across the viscosity transition. In order to trap TOI-216b, α DZ = 10 − 3 in the dead zone requires α MRI (cid:38) 5 × 10 − 2 in the MRI-active zone. If α DZ = 5 × 10 − 4 , α MRI (cid:38) 7 . 5 × 10 − 2 is needed. Conclusions. We describe a new regime of a migration trap relevant for massive (sub-)Neptunes that puts valuable constraints on the levels of turbulent stress in the inner part of their natal disks.


Introduction
Statistical trends in the population of exoplanets tell us that the most frequent planets on close-in orbits fall into the category of super-Earths and sub-Neptunes (with their radii between 1 and 4 R ⊕ ; e.g.Mayor et al. 2011;Howard et al. 2012;Petigura et al. 2013;Winn & Fabrycky 2015).The occurrence rate of these planets increases sharply at around ∼0.1 au (orbital periods ∼10 d; Petigura et al. 2018), which is also where the majority of inner planets in multi-planetary systems reside (Mulders et al. 2018).
Due to the large amount of solids confined within super-Earths (compared to solar system terrestrial planets) and due to the existence of H/He envelopes of sub-Neptunes (e.g.Hadden & Lithwick 2014;Marcy et al. 2014), it is believed that their formation was likely finished before their natal disk was dispersed (e.g.Lambrechts et al. 2019;Ogihara & Hori 2020;Ogihara et al. 2020;Venturini et al. 2020a,b).If so, these planets must have been subject to disk-driven planetary migration which is Movies associated to Figs. 3,10,and 13 are only available at https://www.aanda.orgoften inward-directed (Goldreich & Tremaine 1979;Korycansky & Pollack 1993;Ward 1997;Miyoshi et al. 1999;Masset 2002;Tanaka et al. 2002) and thus an explanation is needed on how to stop the orbital migration and how to explain the aforementioned clustering at 0.1 au.One possibility is the migration trap mechanism of Masset et al. (2006).When a migrating planet encounters a surface density drop, the component of the co-rotation torque driven by the vortensity gradient (i.e. the horseshoe drag; Ward 1991) exhibits a positive boost that can counteract the negative Lindblad torque of the spiral arms.Consequently, the planetary orbit converges to a radius where the total torque is zero and the migration effectively ceases.The clustering at 0.1 au might then be associated with the inner disk rim 1 (e.g.Terquem & Papaloizou 2007;Izidoro et al. 2017;Brasser et al. 2018) where a surface density drop is expected, perhaps due to a viscosity transition from the outer dead zone to the inner zone with A&A 666, A63 (2022) active magneto-rotational instability (MRI; Flock et al. 2017).Flock et al. (2019) present a detailed physical model of the inner disk rim for a typical T Tauri star and predict a migration trap positioned near 0.1 au for a range of model parameters.
A specific example of an observed exoplanetary system that requires the trapping of the inner planet is TOI-216 (Dawson et al. 2019(Dawson et al. , 2021;;Kipping et al. 2019), which was studied in our recent work Nesvorný et al. (2022).The system is comprised of a sub-Solar K-type host star (M = 0.77 M ), an inner Neptuneclass TOI-216b (M p = 0.059 M Jup , P orb 17.1 d), and an outer half-Jupiter TOI-216c (M p = 0.56 M Jup , P orb 34.6 d), with the planet pair captured firmly in the 2:1 mean-motion resonance.By examining various formation scenarios, Nesvorný et al. (2022) argue that trapping TOI-216b close to its observed orbital distance is an essential prerequisite for a successful creation of the 2:1 resonant lock and its subsequent (hydro)dynamical evolution towards the observed state (namely the eccentricity of TOI-216b and the libration amplitude in the resonance).Since the model of Flock et al. (2019) can be applied to TOI-216 without significant changes, Nesvorný et al. (2022) conclude that the planet trap was likely facilitated by the viscosity transition at the inner disk rim.
However, we also point out in Nesvorný et al. (2022) that TOI-216b must have been capable of opening a relatively deep gap and thus it is unclear whether it could have been trapped at the inner disk rim efficiently.The reason for this doubt is that if a planet opens a gap, its co-rotation region becomes emptied and thus the contribution of the positive co-rotation torque is diminished.A reliable answer can only be provided by means of detailed hydrodynamic simulations which were not the main objective of Nesvorný et al. (2022) and we thus focus on them here.
Although we chose TOI-216b as our reference case, the efficiency of the trapping mechanism is a question that applies to sub-Neptunes in general.To demonstrate that, let us recall that the ability of a planet to open a gap and thus switch from the Type I to Type II migration regime predominantly depends on three quantities: the planet-to-star mass ratio q, the disk alpha viscosity α, and the disk aspect ratio h.One of the available predictions for the gap depth (Kanagawa et al. 2018) leads to where Σ gap is the minimal surface density at the gap bottom, Σ 0 is the unperturbed surface density, and The key fact to take into account is that h is very small at 0.1 au.Flock et al. (2019) find h 0.023 near the viscosity transition and this value is relatively well constrained because it is simply a result of the thermal balance in a stellar-irradiated disk2 .Before reaching the inner rim, the planet is embedded in the dead zone where α = 10 −4 -10 −3 can be expected, with the lower limit corresponding, for example, to a hydrodynamic turbulence (Nelson et al. 2013;Klahr & Hubbard 2014) and the upper limit arising from magneto-hydrodynamic (MHD) models (Flock et al. 2017).
Choosing the planet-to-star mass ratio q = 3 × 10 −5 (corresponding to a ten-Earth-mass planet orbiting a Solar-mass star), one obtains 3 Σ gap /Σ 0 0.018-0.15.Such an excavation of the corotation region would affect the co-rotation torque without a doubt (Kanagawa et al. 2018).However, it is not a priori clear if Eqs. ( 1) and (2) still apply once the planet reaches the viscosity transition as they were derived for constant α.Our simulations are thus meant to clarify the picture.
To date, there have not been a large number of hydrodynamic studies dealing with the migration of planets that open moderately deep gaps across the viscosity transition of the inner rim.In Ataiee & Kley (2021), the planet was either too massive (i.e.having a Jupiter mass) or the aspect ratio was too large (typically h = 0.05; the value of h = 0.03 was only considered in a limited number of cases).Romanova et al. (2019) performed detailed 3D simulations at the disk-cavity boundary with the main parameters h = 0.03, q = 1.5 × 10 −5 or 4.5 × 10 −5 , and zero viscosity.They find very efficient trapping, but their simulation timescales were relatively short and the planet was allowed to migrate even before the (possible) gap was relaxed (see their Fig.B2).Moreover, they did not consider a viscosity transition.Finally, Faure & Nelson (2016) studied the planet migration at the inner edge of a dead zone in the presence of a cyclically evolving vortex.For their setup, they identify q ∈ (3 × 10 −5 , 10 −4 ) as the optimal range for efficient trapping.Based on the aforementioned works, it is not straightforward to assess whether (sub-)Neptunes are subject to efficient trapping at the disk edge, and re-investigating the migration trap in this planet mass range is thus worthwhile.
At the same time, the trap should not be too efficient when q increases towards super-Neptunes because the exoplanet occurrence rate steeply decreases already around the Neptune mass and results in the Neptunian desert.Having an efficient migration trap in this mass range might result in an inconsistency between the planet migration theory and observations.We think that this is an important dynamical aspect and it is further elaborated on throughout the paper.We point out, however, that our study cannot rule out other mechanisms that might have created the Neptunian desert, such as photoevaporation (e.g.Owen & Wu 2013;Lopez & Fortney 2013;Owen & Lai 2018;McDonald et al. 2019;Hallatt & Lee 2022), tidal stripping (e.g.Matsakos & Königl 2016;Owen & Lai 2018), or pebble isolation (e.g.Bitsch 2019;Venturini et al. 2020b).
The paper is structured as follows.Our numerical models are described in Sect. 2. Section 3 is devoted to 2D locally isothermal simulations with the viscosity transition at 1 au which enable us to explore a larger portion of the parametric space and get a basic grasp of the processes that play a decisive role for the efficiency of the migration trap.In Sect.4, we perform more realistic non-isothermal simulations at 0.1 au to finalize our findings.A general discussion is provided in Sect. 5 and Sect.6 summarizes our conclusions.

Methods
We used the Fargo3D code (Benítez-Llambay & Masset 2016), which is a finite-difference solver based on upwind methods (Stone & Norman 1992), paired with the fast orbital advection algorithm (Masset 2000).We accounted for a viscous protoplanetary disk (its gaseous component), a central star, and a single planet embedded in the disk.The system was evolved on a 2D polar staggered mesh characterized by the radius r and azimuth θ.The frame was centred on the host star, but co-rotated with the planet.

Locally isothermal 2D model
Governing equations of the model read as follows: where Σ is the gas surface density, u is the velocity vector, P is the pressure, T is the viscous stress tensor, and φ is the gravitational potential.The equation of state is given by where c s is the sound speed, Ω K is the Keplerian frequency, and H is the pressure scale height.The latter is a radially fixed function of the aspect ratio h = H/r, which was uniform in our locally isothermal model (i.e. the disk flaring was neglected).
The potential component related to the planet was smoothed via where d is the distance from a grid cell to the planet and r sm = 0.6H is the smoothing length that allows one to avoid numerical divergence at d = 0 and also leads to disk torques that resemble those found in 3D simulations (Müller et al. 2012).The gravitational potential also included disk-and planet-driven indirect terms that arise from the usage of a non-inertial reference frame.
Since we are mostly concerned with planets that open (at least partial) gaps, it is likely that these planets also accumulate gas in a circumplanetary disk.However, this region is underresolved and due to the lack of self-gravity, its contribution to the torque exerted on the planet is incorrect.Therefore, it is customary to exclude a part of the gas contained within the Hill sphere when calculating the disk torque (but not when computing the back-reaction of the planet on the gas).When evaluating the gravitational acceleration exerted by a gas cell on the planet, we multiplied by the tapering formula adopted from Crida et al. (2008): where p = 0.8.In addition, we enabled the module of Fargo3D that activates the correction for the shift of Lindblad resonances that is also related to the lack of self-gravity (Baruteau & Masset 2008).The location of resonances is improved by subtracting the azimuthally averaged gas density from Σ prior to the torque evaluation.
The disk was initialized with a uniform mass accretion rate where the kinematic viscosity ν was calculated using the alpha prescription (Shakura & Sunyaev 1973) To mimic the viscosity transition and the surface density drop at the inner disk rim, we used (for a similar temperature-dependent transition, see Eq. ( 26) or Flock et al. 2016) where α MRI accounts for the alpha viscosity in the MRI-active zone inwards from the transition radius r MRI and α DZ characterizes the viscous stress in the dead zone at r > r MRI .The sharpness of the viscosity transition was controlled by the parameter ∆r MRI .The radial velocity field was initially and the azimuthal velocities v θ were set by the centrifugal balance between the gravity of the star and the pressure support within the disk.The initial conditions for Σ, v r and v θ were also used as boundary conditions.To provide a smooth transition between the boundaries and the perturbed disk, we damped Σ and v r towards their initial values within radially narrow zones adjacent to the boundaries.The damping recipe was that of de Val-Borro et al. (2006).

Non-isothermal 2D model
We relaxed the locally isothermal approximation from previous Sect.2.1 by replacing Eq. ( 5) with where γ is the adiabatic index and E is the internal energy density.To describe the temporal evolution of E, we followed the one-temperature approximation (i.e. a thermodynamic equilibrium between the gas and photons of thermal radiation) and considered (e.g.Kley & Crida 2008;Pierens 2015;Chrenko et al. 2017;Ziampras et al. 2020) where Q visc is the viscous heating term, Q irr is the heating from stellar irradiation, Q vert is the vertical cooling by radiation escape from the disk, and the term −P∇ • u accounts for the compressional heating.Let us note that although our implementation allows for in-plane radiative diffusion, it was not considered in this work and therefore omitted in Eq. ( 13).
The closure relation for viscous heating is given in Appendix A (Eq. A.4).The irradiation-cooling balance can be written as where σ is the Stefan-Boltzmann constant, T irr is the midplane temperature of a passive disk, and T is the instantaneous temperature related to the internal energy density as E = ρc V T , c V being the heat capacity at constant volume.To calculate the effective optical depth τ eff which allowed us to estimate local radiative losses, we followed Hubeny (1990) andD'Angelo &Marzari (2012): with In writing Eq. ( 16), we introduced a single grey opacity κ, which was considered uniform.
A63, page 3 of 20 A&A 666, A63 (2022) Ultimately, we aimed for our model to be applicable at the inner rim of protoplanetary disks, which is thermodynamically complex.Flock et al. (2016Flock et al. ( , 2017) ) performed vertically resolved radiation (magneto)hydrodynamic studies of the inner rim structure and identified four distinct regions when moving from the star towards larger radii (see also Fig. 11): (A) an optically thin region of pure gas, (B) an optically thin halo with a low dust concentration, (C) a rounded off condensation front of dust grains that is optically thick to stellar irradiation (and can be followed by a radially narrow self-shadowed region; Flock et al. 2019), and (D) a standard flaring disk with the photospheric height of several H.The viscosity transition at the inner rim is related to the ionization temperature that triggers the MRI and it is often located in the proximity of the boundary between regions C and D. Therefore, our goal was to reproduce the thermal structure of regions C and D4 because they contain the first possible trap5 that an inward-migrating planet encounters when spiraling towards the star.
Although our model is not vertically resolved, it is possible to follow the approach of Ueda et al. (2017), who analytically studied the thermal balance found in numerical simulations of Flock et al. (2016), and recover the temperature profile in regions C and D with a satisfactory precision (as we show in Sect. 4 and Fig. 11).We took where T trans is the temperature in region C and T thick is the temperature in region D. For the former, Ueda et al. (2017) find where r cond is the radius of dust condensation at all heights above the midplane and T cond is the temperature of the front of fully condensed dust.The remaining closure relations are where and In these expressions, r rim is the innermost radius where the disk midplane becomes thick to stellar irradiation, is the emission-to-absorption efficiency of dust grains, T is the stellar temperature, T ev is the temperature in the dust halo where the grains are evaporating (region B), R is the stellar radius, f ph,cond is the photospheric height at r cond (given in multiples of H), and M is the stellar mass.
For the temperature in region D, we derive a slightly modified formula with respect to Ueda et al. (2017) in Appendix B: where f ph is the photospheric height in region D (again given in multiples of H).
In order to mimic the decrease of the vertical optical depth when moving from region D to region C, we introduced an ad hoc modification of τ opt by multiplying it with where f τ,trans = 0.1 and f τ,thick = 1.Equipped with all the formulae for the heating and cooling terms, Eq. ( 13) was solved for temperature T in an implicit form using E = ρc V T .As we omitted the in-plane radiative diffusion, the implicit update has a simple analytic solution that can be performed cell by cell.
In the non-isothermal model, we used the adiabatic versions of the sound speed and pressure scale height: All planet-related calculations and the smoothed planetary potential remained the same as in Sect.2.1, the only difference was that we used the azimuthally averaged H when computing the potential smoothing length.Equations (8-11) remain valid as well.But finding the initial state of the disk required a more detailed approach because ν depends on the temperature via c s and this dependence propagates into the calculation of Σ from Ṁ.We thus proceeded iteratively: Starting from an initial guess of the radial temperature profile, we calculated the density profile from Eqs. (8-10).Then we kept the density fixed and we evolved Eq. ( 13) over a time interval ∆t that was similar to the characteristic radiation diffusion timescale in the disk.The steps were repeated until the relative change of T and Σ became smaller than 10 −5 .Subsequently, we confirmed that a correct equilibrium state was found by propagating the full set of Eqs.(3), (4) and (13) over 1000 orbital periods P orb (measured at the viscous transition).During the relaxation, the disk was effectively kept in 1D by assuming an axial symmetry and setting the number of azimuthal sectors equal to one.After the relaxation, we prepared the disk for simulations with an embedded planet by expanding all hydrodynamic fields in azimuth and by converting velocities to the frame co-rotating with the planet.
The boundary conditions were constructed as follows.First, we used a power-law extrapolation of T from the active grid cells to the ghost cells6 .Next, we used T in the ghost cells to find ν.Afterwards, v r was set in the ghost cells using Eq. ( 11).The calculation of Σ differed for the inner and outer boundary in radius (see also Bitsch et al. 2014).At the inner boundary, we kept the inward-directed mass flux through the boundary uniform by using Σ g = Σ a ν a /ν g , where the subscripts 'g' and 'a' stand for the ghost and active cells, respectively.At the outer boundary, we used Σ g = Ṁ/(3πν g ) in order to replenish the disk material in a uniform way.Finally, E in the ghost cells was recovered from T and Σ.
Damping zones were used in the non-isothermal model as well.During the relaxation stage, we damped only v r so that Σ could adjust and reach an equilibrium.The target value of v r , towards which we damped, was again given by Eq. ( 11), which depends on the instantaneous local temperature.During the simulations with an embedded planet, we damped Σ and v r towards the equilibrium disk state found at the end of the relaxation stage.

Locally isothermal simulations at 1 au
In the first part of our study, we perform 2D locally isothermal simulations (Sect.2.1) with the viscosity transition positioned roughly at r MRI = 1 au.Locating the viscosity transition further out (with respect to 0.1 au) means that a simulation covering tens of thousands of orbital timescales translates to physical times that are comparable to the expected migration timescales in the Type II regime.In other words, we are able to cover the migration history of gap-opening planets as they approach the viscosity transition and as they interact with it.Using the locally isothermal approximation, and thus neglecting the disk thermodynamics, is not very realistic yet it simplifies the parametric space and it again helps to reduce numerical demands.In the end, the simplicity can be advantageous because it can provide hints on the robustness of our results.
Table 1 defines parameters that we vary in our simulations.They are specified throughout the following sections.Table 2, on the other hand, summarizes parameters that are kept fixed.As we intend to use TOI-216b as our reference case, we adopt the mass of its host star M = 0.77 M and, unless specified differently, the planet-to-star mass ratio is q = 7.3 × 10 −5 .The value of h = 0.023 is motivated by the aspect ratio at the inner disk rim (see also Fig. 11).As for the computational domain, it is chosen such that the Hill radius of TOI-216b is resolved by 7 cells and the half-width of the horseshoe region (e.g.Jiménez & Masset 2017) by 19 cells.The radial spacing is logarithmic and the grid resolution is such that individual cells remain roughly square-shaped to reduce numerical anisotropy.

The role of the α-viscosity contrast for trapping TOI-216b
We start by placing TOI-216b at a p (t 0 ) = 1.6 au, outwards from the viscosity transition.The viscosity transition width parameter is set to ∆r MRI = 2H = 0.046 au.Two values are explored for the dead-zone viscosity, α DZ = 5 × 10 −4 and 10 −3 .For each of these α DZ values, we vary the MRI-active viscosity as α MRI = (10 −2 , 2.5 × 10 −2 , 5 × 10 −2 , 7.5 × 10 −2 , 10 −1 ).Since the disk-driven migration scales with the disk mass, it would be ideal to maintain the unperturbed Σ(r) profile the same in all our simulations.However, this can never be fulfilled over the whole extent of the disk when varying the viscosities and initializing the disk via Eq.( 8) at the same time.We decided to maintain the same Σ at least in the dead zone of the disk, which is the region that (i) controls the approach of the planet towards the viscosity transition and (ii) contains most of the disk mass in the given radial range.This is done by setting Ṁ = 10 −8 M yr −1 when α DZ = 10 −3 and Ṁ = 5 × 10 −9 M yr −1 when α DZ = 5 × 10 −4 .The planet mass is gradually increased over the first 1000 orbital timescales and the planet itself is held in a fixed circular orbit for 10 000 and 15 000 orbital timescales when α DZ = 10 −3 and 5 × 10 −4 , respectively.Only after a gap was formed and its opening slowed down did we let the planet migrate for 50 000 P orb .Since a p (t 0 ) = 1.6 au, our setup allows the gap to be first opened in the initially smooth part of the dead zone, without any influence of the viscosity (and surface density) transition.This way, the planet first acquires a correct migration rate before it approaches the viscosity transition where the gap profile is expected to become modified.We recall here that a gap-opening planet needs to be allowed to migrate for some time before its migration rate converges (Dürmann & Kley 2015;Robert et al. 2018;Scardoni et al. 2020;Kanagawa & Tanaka 2020;Lega et al. 2021).The reason is that the gap of a non-migrating planet represents a different barrier (although leaky) for the background gas flow than a gap of a migrating planet (Lubow & D'Angelo 2006;Duffell et al. 2014).
Temporal evolution of the semi-major axis a p after the release of TOI-261b is shown in Coloured curves correspond to different simulation times t as indicated by the plot legend.The planet is kept fixed for t = 10 4 and 1.5 × 10 4 P orb when α DZ = 10 −3 and 5 × 10 −4 , respectively.At larger times, the planet is allowed to migrate freely.The vertical dashed line marks the final location of the planet in simulations where trapping occurs.We point out that by setting Ṁ based on α DZ (see Sect. 3.1), the unperturbed density profile in the dead zone is identical in each simulation.
transitions.At first, the migration proceeds as expected for a gap-opening planet.For a short time, the planet experiences inward migration uncoupled from the viscous flow in the disk due to the initial inbalance of the torques (as reflected by the slightly steeper slope at the very beginning of the migration curves).But then the migration curves level at a slope tied to the viscous flow (similarly to e.g.Robert et al. 2018).
When the planet migrates across a p = 1.2-1.3au, the part of the inner spiral arm that is propagating through the decreased density region below the viscosity transition becomes progressively larger.Consequently, the outer spiral arm takes over and the Lindblad torque becomes more negative, resulting in an episode of fast inward migration in the vicinity of the viscosity transition.What happens next depends on the viscosity values used to model the transition.When α DZ = 10 −3 , the planet is trapped at the viscosity transition for α MRI ≥ 5 × 10 −2 , otherwise it continues migrating inwards.When α DZ = 5 × 10 −4 , the trapping is only possible for α MRI ≥ 7.5 × 10 −2 .Such a dependence on the α-viscosity suggests an interplay between the gap depth and the trap efficiency -as the gap carved by the planet gets deeper in less viscous dead zones, larger α MRI is then required to (partially) refill the gap so that the trapping mechanism can still work.Therefore, we examine the gap profile next.

Evolution of the gap profile
Figure 2 shows the temporal evolution of the radial surface density profile Σ(r) for four selected simulations that combine the minimum and maximum α-viscosities from our parametric range.To obtain the Σ(r) profiles, we calculated the arithmetic mean from Σ(r, θ) in each radial ring of the domain.We also immediately see a substantial difference in the gap profile between simulations in which the planet is trapped and those in which it continues migrating inwards.For simulations where the trapping fails (bottom row in Fig. 2), we see that the gap mostly retains its characteristic shape even when the planet is crossing the viscosity transition -the gap centre still appears to be close to an inner-outer symmetry, with a small central peak that reflects the mass accumulation in the Hill sphere of the planet.The gap walls do exhibit an inner-outer asymmetry, with the outer wall becoming steeper.As a result, the outer spiral arm contains more gas mass relative to the inner spiral arm and thus drives a strong negative Lindblad torque that is responsible for the episodes of fast inward migration in Fig. 1.When the trapping occurs (top row in Fig. 2), on the other hand, the gap profile becomes modified.The inner half of the gap cannot be distinguished from the unperturbed density profile and a local density peak occurs at the planet location.The outer half of the gap remains cleared of gas, although it is shallower compared to simulations in which there is no trapping.Due to the drop in the surface density across the viscous transition, the outer spiral arm inevitably dominates over the inner one in this case as well.We thus need to focus our attention on the changes in the co-rotation region in order to understand why the torque felt by the planet becomes zero.To do so, the next section analyses 2D maps of the surface density as well as gas flow properties.Fig. 4. As in Fig. 1, but for fixed α DZ = 10 −3 and different values of the planet-to-star mass ratio q (see the plot legend).Two cases of the MRI-active viscosity are considered, namely α MRI = 5 × 10 −2 (top) and 2.5 × 10 −2 (bottom).Dashed lines correspond to simulations in which the torque is evaluated from the whole disk, without any exclusion of the Hill sphere material.Solid lines are affected by the Hill cut (see Eq. ( 7)).For comparison purposes, the case of TOI-216b is labelled and marked with an arrow.

Surface density and gas flow evolution in the co-rotation region
viscosity transition.The gap is regular as well as the streamlines within.However, as the planet continues approaching the viscosity transition, the gas density distribution at the inner gap edge turns into a radially narrow bump.This can be seen near 1.1 au in panel b where the initial step-like profile of Σ(r) turns into a bump that is bordered by the low-density high-viscosity region on the inside and by the planet-induced gap on the outside.The density bump becomes prone to the Rossby wave instability (RWI; Lovelace et al. 1999;Li et al. 2000) that manifests itself via a high-m mode and excites multiple vortices.The instability does not persist as the gap follows the inwardmigrating planet and reduces the gas accumulation at 1.1 au even more, as shown in panel c.Nevertheless, a less prominent local density bump is still maintained within the co-rotation region of the planet.The streamlines in the co-rotation region develop two librating islands, one positioned in the front horseshoe region (leading the orbital motion of the planet) and one located in the rear horseshoe region (trailing the orbital motion of the planet).We notice a front-rear asymmetry, with the front island being azimuthally narrower.Over time, both islands of librating streamlines accumulate gas but the front one is more efficient as indicated by panel d.
Eventually, the co-rotation streamlines reconfigure when the rear overdensity propagates through the co-rotation region in a retrograde sense (with respect to the planet) and merges with the front overdensity.The state at the time of the merger is depicted in panel e.At this t mig , the migration trap is first detected (da p /dt becomes non-negative for the first time).The overdensity in the front horseshoe region becomes a Rossby vortex and launches its own set of spiral arms.The final state of the simulation is shown in panel f.The vortex in the front libration island persists and its peak density (with respect to panel e) increases.
Overall, it is clear that the trapping is facilitated by the co-rotation torque but its nature is different compared to Masset et al. (2006).According to Masset et al. (2006), planets stop because the strong vortensity gradient at the surface density transition modifies the horseshoe drag (Ward 1991).Here, instead, the partial gap opening reduces the initial gas density and changes its shape from a step-like to a bump-like.Thus the vortensity gradient is diminished as well as the actual amount of gas in the co-rotation region that can exert the torque.The trapping is only possible due to the front-rear asymmetry that occurs in the co-rotation region.The front overdensity dominates and exerts a positive torque on the planet (as any azimuthal overdensity leading the orbital motion of the planet would) that balances the negative torque of spiral arms.Consequently, we find that the migration stops.
Although the nature of the trap for TOI-216b is now clear, it is not a priori obvious what is the main driving mechanism for the occurrence of islands of librating streamlines.On one hand, we detect vortices in panels b, e and f of Fig. 3 and those would suggest that the RWI dominates.On the other hand, panels c and d already exhibit the libration islands but there is no apparent sign of vortices and that would suggest that the islands are a natural feature of a partially emptied co-rotation region that overlaps a viscosity transition.To gain further insight concerning this problem, we explore the behaviour of the planet trap for various planet-to-star mass ratios q in the upcoming section.

Varying the planetary mass
In previous sections, the planetary mass was fixed to that of TOI-216b.In this section, we investigate the planet-to-star mass ratios 7 q = (1.5, 3, 4.5, 6, 9, 18, 36) × 10 −5 .Fig. 4 shows the temporal evolution of the planetary semi-major axes a p (t) for fixed α DZ = 10 −3 and two values of the MRI-active viscosity α MRI = 2.5 × 10 −2 and α MRI = 5 × 10 −2 .We selected these two viscosity values because they represent the divide between trapping and not trapping TOI-216b (see the labelled green curve in Fig. 4).The remaining parameters remain the same as in Sect.3.1 apart from the total simulation timescale after the release of the planet, which we reduced to 40 000 P orb .In terms of the trapping efficiency, Fig. 4 reveals that all studied planets less massive than TOI-216b become trapped at the viscosity transition while more massive planets avoid the trap and continue migrating inwards (with one peculiar exception mentioned below).
It is also worth noting that planets with q ≤ 4.5 × 10 −5 experience a runaway migration over the time interval 100-350 P orb that is barely resolved in Fig. 4. To exclude the possibility that the runaway migration occurs due to the exclusion of some of the gas in the Hill sphere from the torque evaluation, which might not be fully justified for low values of q, we recalculated the two least massive cases while considering the torque exerted by the whole disk (dashed lines in Fig. 4).The Hill cut has a minimal influence on the semi-major axis evolution.The reason for the fast inward migration is that these planets open partial gaps and then they enter the Type III migration regime (Masset & Papaloizou 2003;Pepliński et al. 2008a).To our knowledge, this is the first numerical suggestion that super-Earths or mini-Neptunes might undergo inward Type III migration under the conditions specific for the inner disk region (mainly the very low aspect ratio h and a suitably large Σ).Despite their runaway migration, these 7 Our range of q would translate to planetary masses M p = (3.85, 7.7, 11.55, 15.4, 23.1, 46.2, 92.4) M ⊕ in the TOI-216 system or to M p = (5,10,15,20,30,60,120)  ] radial distance r [au] q = 9 x 10 -5 t/P orb = 9000 Fig. 6.Radial surface density profiles for cases discussed in Fig. 5 (planet-to-star mass ratios q are given by labels in individual panels).We plot the initial azimuthal average (solid grey curve), the azimuthal average at t orb = 40 000 P orb for trapped planets and at t orb = 9000 P orb for q = 9 × 10 −5 (black solid curve), azimuthal averages taken over the span of the front (θ > 0) and rear (θ < 0) horseshoe regions (red and blue solid curves, respectively), and radial cuts Σ(r, θ = π/3) and Σ(r, θ = −π/3) (orange and purple dashed curves, respectively).We point out that the Hill sphere material is excluded from calculations of all azimuthal averages.The instantaneous position of the planet is marked by the black arrow.Unlike in Fig. 2, the vertical scale is linear.The spread of the curves (notably the orange and red ones) indicates the development of the front-rear asymmetry in the co-rotation region.front libration island; (iv) a front-rear asymmetry in the gas density of the co-rotation region develops; (v) spiral arms launched by a front vortex are apparent for q ≥ 6 × 10 −5 .
In Fig. 5, planets with q ≤ 7.3 × 10 −5 are trapped.Based on this fact and based on the description provided in the previous paragraph, we conclude that the discussed simulation set captures a transition from the classical trap of Masset et al. (2006) (that should dominate e.g. for q = 1.5 × 10 −5 for which we do not see any strong front-rear asymmetry in the gas distribution) to a trap driven by the front-rear asymmetry of the co-rotation region for planets with moderate gaps (e.g.q = 4.5 × 10 −5 -7.3 × 10 −5 ) and finally to untrappable planets that clear their co-rotation regions too efficiently (e.g.q = 9 × 10 −5 ).
To further support these claims, Fig. 6 shows the radial surface density profiles for the cases discussed in Fig. 5.We can see in detail how the initial step-like transition of Σ(r) changes to a bump-dominated profile, whose amplitude decreases with increasing q.Additionally, by plotting the azimuthal averages over front and rear horseshoe regions and also radial cuts through the domain at θ ± π/3, we reveal the development and magnitude of the front-rear asymmetry of the horseshoe region.The most prominent asymmetry occurs in the range q = 4.5-7.3× 10 −5 and for these planets, it produces the positive torque necessary for their trapping.For q = 9 × 10 −5 , the failure of the migration trap is caused by the decreased level of the front-rear asymmetry and also by the overall emptying of the co-rotation region.
Finally, let us return to the question of the relative importance of the RWI compared to planet-induced perturbations of the gas flow (see the end of Sect.3.3).To proceed with the analysis, Fig. 7 shows maps of the inverse potential vorticity (or inverse vortensity) defined as where the vorticity component is perpendicular to the disk plane and the term related to the orbital frequency translates ζ −1 to the inertial frame.The inverse potential vorticity has several diagnostic applications in our case.First, we look for the local maxima of ζ −1 to trace vortices that are the usual outcomes of non-linear RWI growth8 .Second, the vortensity-driven horseshoe drag scales as Γ HS ∝ Σd log ζ −1 /d log r (Masset et al. 2006) so that steep radial gradients of ζ −1 can help us judge if the classical migration trap of Masset et al. (2006) still applies or not.
In Fig. 7, we focus on three cases q = (1.5, 4.5, 7.3) × 10 −5 .For the lowest planet mass, we see that a steep positive gradient of ζ −1 across the co-rotation is maintained and that confirms that the planet is trapped primarily due to the mechanism of Masset et al. (2006).Additionally, we notice a small asymmetry at the outer separatrix of the horseshoe region -its front part has a slightly larger ζ −1 compared to the rear part.This is due to two effects: (i) the vorticity production across the shock front of the outer spiral arm (Koller et al. 2003) that propagates through a low-viscosity dead zone where the ability of the planet to shock the disk is higher (Rafikov 2002) and (ii) the low-density material being carried from the inner disk towards the outer disk by rear horseshoes (Romanova et al. 2019).Both effects (growing ∇ × u and mixing with a lower Σ) reduce ζ −1 along the outer separatrix of the rear horseshoe region.
The middle and bottom panels of Fig. 7 show that the gradient of ζ −1 across the co-rotation is reduced with increasing planet mass, as we already anticipated.Moreover, we notice that a steep isolated local maximum of ζ −1 only appears for q = 7.3 × 10 −5 .For q = 4.5 × 10 −5 , the peak values do not exceed those found  30)).The colour scale is saturated to highlight features in the co-rotation region.
for q = 1.5 × 10 −5 and thus the planet with q = 4.5 × 10 −5 does not develop a vortex in its co-rotation region (which is consistent with the lack of vortex-driven spiral arms for this planet mass).Yet the co-rotation region of q = 4.5 × 10 −5 already exhibits a front-rear asymmetry dominated by the front island of librating streamlines.In addition, we recall that the azimuthally averaged density bump across the co-rotation is actually larger for q = 4.5 × 10 −5 than for q = 7.3 × 10 −5 (panels c and e, respectively, We point out that the total width of the viscosity transition (as marked by the shaded area) is not exactly equal to ∆r MRI ; the influence of the parameter is described by Eq. ( 10).
in Fig. 6).The fact that the sharper peak in gas density remains Rossby-stable means that the vortex found for q = 7.3 × 10 −5 has to be a secondary effect launched by the streamline perturbations due to the influence of the planet.This answers the question on the relative importance of these two effects: the islands of librating streamlines develop regardless of the RWI but the RWI can be triggered secondarily if an island accumulates enough mass or exhibits a substantial change in the curl of the gas velocity.Before moving on, let us point out that although we studied the dependence on q in this section, the described effects are related rather to the gap-opening ability of a planet that scales according to Eqs. ( 1) and (2).Therefore, similar effects can be achieved by variations of α DZ and α MRI (which is consistent with results reported in Sect.3.1) or by varying h (not explored in this paper).

Varying the viscosity transition width
In this section we modify the parameter ∆r MRI that controls the width of the viscosity transition.We set α DZ = 10 −3 , α MRI = (10 −2 , 2.5 × 10 −2 , 5 × 10 −2 , 7.5 × 10 −2 , 10 −1 ) and we explore two cases with ∆r MRI = H = 0.023 au and ∆r MRI = 4H = 0.092 au.The remaining parameters follow Sect.3.1 with the exception of the total simulation timescale after the release of the planet, which we reduced to 30 000 P orb .
The migration tracks that we obtained are shown in Fig. 8.We see that if the viscosity transition is sudden, the range of viscosities in the MRI-active zone that facilitate the trapping increases to α MRI ≥ 2.5 × 10 −2 compared to Our exploration of the parametric space does not go beyond the simple experiment shown in Fig. 8 yet it reveals that the width of the viscosity transition might be equally important for the efficiency of the migration trap as the actual values of the viscosity itself.For a reference, we recall that Flock et al. (2017) found rather sharp viscosity transitions occurring over 1 or 2H but their Ohmic resistivity model was simplified.Future magnetohydrodynamic studies should therefore improve the ionization profile and disk chemistry to carefully reassess the radial disk range over which the MRI is triggered.

Exploring the static torque
In previous sections, we saw that the migration trap for sub-Neptunes and Neptunes at the inner rim is caused by the frontrear asymmetry of the co-rotation region.Similar asymmetries, although of a different origin, often occur as a result of planet migration that leads to streamline deformation and, since they would not occur for planets in fixed orbits, the resulting torques are usually referred to as dynamical torques.Here, we perform simulations of TOI-216b in a fixed circular orbit and measure the static torque (i.e. the torque acting on a non-migrating planet) for a set of semi-major axes a p distributed in a close proximity of the viscosity transition.Our aim is to check whether the static torque measurements can also recover the trap and its position compared to simulations with a migrating planet (Sect.3.1).In other words, we want to check whether the trapping mechanism depends on the migration history of the planet or not, which is an important question.
The parameters remain the same as in Sect.3.1 but we limit the explored viscosities to α DZ = 10 −3 and α MRI = (2.5 × 10 −2 , 5 × 10 −2 , 10 −1 ).Additionally, the planetary semi-major axis remains fixed and we take a p distributed evenly between 0.975 au and 1.15 au with a step of ∆a p = 0.025 au.The planetary mass is again increased over the first 1000 P orb and each simulation covers 20 000 P orb to allow the disk torque to converge.Figure 9 shows our measurements of the static torque Γ normalized to Γ 0 = Σ p a 4 p (Ω p q/h p ) 2 (e.g.Korycansky & Pollack 1993), where the subscript refers to the planet location and the surface density is that of an unperturbed disk.Wherever Γ/Γ 0 > 0 in Fig. 9, outward migration is expected and vice versa.For α MRI = 10 −1 , linear interpolation between the static torque measurements predicts trapping at a p = 1.1095 au compared to 1.1009 au found for the migrating planet.For α MRI = 5 × 10 −2 , the zero-torque radius is located at a p = 1.0834 au compared to the dynamical trap location 1.0708 au.For α MRI = 2.5 × 10 −2 , the static torque is always negative and the planet would not become trapped.
Overall, Fig. 9 is in a very good agreement with Fig. 1.The trapping is predicted for the same values of α MRI and the trap location differs only by 0.01 au on average.Although we advocated in Sect.3.1 that the planet needs to be migrating for the Type II torque to be evaluated correctly, we demonstrated here that static torque measurements are applicable when searching for the existence of the migration trap and its location.
Finally, we check that the surface density distribution and streamlines behave in a similar manner compared to experiments with a migrating planet, mainly in the terms of the front-rear asymmetry of the co-rotation region.Figure 10 depicts the final state of our simulation with α DZ = 10 −3 , α MRI = 10 −1 and TOI-216b fixed at a p = 1.1 au, close to the expected trap location.By comparing to panel f of Fig. 3, we immediately see that the disk response to planet-driven perturbations is qualitatively the same and the torque felt by the planet is inevitably affected by the overdensity in the front horseshoe region.

Non-isothermal simulations at 0.1 au
The aim of this section is to improve the realism of our simulations and to support or refine our findings from Sect. 3. We thus study a non-isothermal disk with its thermodynamics driven by compressional, viscous and irradiative heating balanced by simplified vertical radiative cooling (Sect.2.2).At the same time, we shift the inner disk rim inwards, to a location close to the current orbit of TOI-216b, r MRI = 0.11 au (Nesvorný et al. 2022).The A63, page 12 of 20  remaining fixed parameters of our non-isothermal simulations are summarized in Table 3.

Model verification
In order to demonstrate the reliability of our 2D non-isothermal model, we compared the equilibrium radial profiles of several characteristic quantities to those found in vertically resolved 3D simulations with radiation transfer.The comparison is shown in Fig. 11.The displayed case is for Ṁ = 10 −9 M yr −1 , α DZ = 10 −3 , α MRI = 5 × 10 −2 .The reference 3D simulation in Fig. 11 is based on the radiative hydrostatic approach of Flock et al. (2016Flock et al. ( , 2019)), which we have recently re-implemented in Fargo3D (Chrenko et al., in prep.).The viscosity transition in the 3D reference simulation is temperature-dependent, following To model the radius-dependent α-transition in the 2D simulation (see Eq. 10), we chose r MRI = 0.11 au together with ∆r MRI = 2.7H 0.007 au and Fig. 11 (top panel) shows that the resulting profile approximates the α(T ) profile of the reference 3D simulation well enough.Consequently, the profile of Σ(r) is also recovered with good accuracy.Let us also point out that the fiducial width of the viscosity transition used in Sect.3, ∆r MRI = 2H, was slightly narrower than we assume here but still relatively similar.
As for the temperature profile (which directly sets the aspect ratio profile), bottom panel of Fig. 11 reveals that we obtain the correct radial dependence in the optically thick flared disk (region D) as well as the steep increase and the plateau near the rounded off evaporation front of silicates (region C).The inner part of the temperature plateau is too extended but this is because we do not attempt to model the dust halo (region B in Fig. 11) with our 2D non-isothermal model (as explained in Sect.2.2).This inconsistency does not influence our results because whether the planet becomes trapped at the viscosity transition or not is decided at larger radii.

Parametric study of the static torque
In our non-isothermal simulations, we focus mainly on static torque measurements motivated by the success of this method in Sect.3.6.Due to computational limitations, we opt for a slightly coarser grid resolution compared to Sect. 3 -for TOI-216b (q = 7.3 × 10 −5 ), we resolve the Hill radius and the half-width of the horseshoe region by 5 and 13 cells, respectively.To explore the free parameters, we combine three values of α DZ = (2.5, 5, 10) × 10 −4 , four values of α MRI = (2.5, 5, 7.5, 10) × 10 −2 and three values of q = (5, 7.3, 9) × 10 −5 .The orbital radius of the planet varies as a p = (0.105, 0.11, 0.115, 0.12, 0.125, 0.13) au.We again scale Ṁ in the same way as α DZ but, compared to locally isothermal simulations, we cannot guarantee an identical surface density profile in the dead zone.The reason is that modifying Ṁ leads to a different Q visc and thus the resulting local temperature becomes slightly different, as well as ν and Σ.However, the differences between individual simulations are rather marginal because Q irr dominates over Q visc for the chosen parameters.We introduce the planet over 500 P orb (r MRI ) and the total timescale of each simulation is 2500 P orb (r MRI ).
Figure 12 summarizes our sweep of the parametric space in terms of the normalized torque γΓ/Γ 0 .The migration trap is expected to occur wherever the torque switches from negative to positive when moving from larger orbital radii inwards.The general trend in Fig. 12 is such that the trapping becomes less efficient with increasing q (left to right in Fig. 12), decreasing α DZ (top to bottom in Fig. 12) and decreasing α MRI (when the torque is positive, it has a larger magnitude e.g. for crosses than for triangular data points).In other words, for parameters that allow opening of a more prominent gap in the dead zone (i.e.decreasing α DZ or increasing q), a stronger viscosity contrast A63, page 13 of 20 A&A 666, A63 (2022) α DZ = 2.5 x 10 -4 Fig. 12. Normalized static torque γΓ/Γ 0 measured in our non-isothermal simulations as a function of the planetary orbital distance a p .We show measurements for q = 5 × 10 −5 (left column), 7.3 × 10 −5 (TOI-216b; middle column) and 9 × 10 −5 (right column).Individual rows correspond to different viscosities α DZ in the dead zone (see the boxed labels) and colour-coded points with line segments correspond to different viscosities α MRI in the MRI-active zone (see the plot legend).The horizontal dashed line is where the total torque is zero.We point out that the extent of the vertical axis in the left column differs from the middle and right columns.
over the transition zone is required for the trap to operate.If the gap is too deep, the trap fails even for the largest tested value of α MRI = 10 −1 .Looking at the case of q = 5 × 10 −5 , which for M = 0.77 M translates to M p = 13 M ⊕ , we see that it can be trapped relatively easily.The only combination of viscosities that fails to trigger the trap is α DZ = 2.5 × 10 −4 and α MRI = 2.5 × 10 −2 .However, based on the general trend described in the previous paragraph, we can predict that it would get progressively difficult to trap the planet if the dead zone viscosity was lower than considered in our parametric sweep, e.g.α DZ 10 −4 .
Focussing now on TOI-216b with q = 7.3 × 10 −5 , we see that α DZ = 10 −3 requires α MRI 5 × 10 −2 in the MRI-active zone to sustain the trap and α DZ = 5 × 10 −4 leads to trapping only if α MRI 7.5 × 10 −2 .For α DZ = 2.5 × 10 −4 , the trap does not exist and the planet would continue migrating inwards.Since trapping of TOI-216b at the viscosity transition is needed to explain the system's architecture (Nesvorný et al. 2022), the results obtained here put constraints on the viscosity values in the central parts of the natal disk of TOI-216.
Finally, we investigate the surface density distribution and streamlines.For a comparison with non-isothermal simulations, we select the viscosity combination of α DZ = 10 −3 and α MRI = 5 × 10 −2 .Figure 13 shows the density maps for planets positioned at a p = 0.12 au.We chose this combination of parameters to (i) provide a straightforward comparison with Fig. 5 and (ii) cover various total torques (from Fig. 12, we see that the torque is large and positive for q = 5 × 10 −5 , small and positive for TOI-216b, small and negative for q = 9 × 10 −5 ).
In general, Fig. 13 is in accordance with Sect. 3. We recover (i) the front island of librating streamlines, (ii) the mass accumulation within this island that can counteract the negative Lindblad torque, and (iii) the overall depletion of the co-rotation region with increasing q.The only difference that we notice is in the general behaviour of streamlines in the co-rotation regiona great part of these streamlines goes from the inner to the outer disk without encountering the planet.Nevertheless, we confirm here that the migration trap at the inner disk rim for planets with moderate gaps is caused by the presence of a prominent overdensity located in the front horseshoe region, in accordance with our locally isothermal simulations.

Implications
Perhaps the most intriguing implication of our study is that sub-Neptunes and Neptunes with q similar to our parametric span can only be trapped at the inner disk rim for specific combinations of viscosities in the inner MRI-active zone and the outer dead zone.Therefore, we obtain hints about processes responsible for generating turbulent stress in natal disks of exoplanetary systems that harbour said sub-Neptunes and Neptunes at 0.1 au.We recall that the trapping is required for both in-situ (e.g.Chatterjee & Tan 2014;Jankovic et al. 2019) and migratory scenarios of planet formation (e.g.Ida & Lin 2010;Cossou et al. 2013;Izidoro et al. 2017) because without the disk torque becoming zero, a concentration of planets would never be achieved (see also Lambrechts et al. 2019;Venturini et al. 2020a,b).
The constraints that we find are still somewhat broad (for instance, trapping of TOI-261b requires α MRI 5 × 10 −2 when α DZ = 10 −3 or α MRI 7.5 × 10 −2 when α DZ = 5 × 10 −4 ) yet they are valuable because (i) they confirm that the MRI activity in the central part of the disk needs to be relatively large and (ii) the "dead zone" cannot become fully dead but needs to maintain a considerable residual level of viscosity (e.g.α DZ ≥ 5 × 10 −4 for A63, page 14 of 20 Fig. 13.Logarithm of the gas surface density in non-isothermal simulations with α DZ = 10 −3 and α MRI = 5 × 10 −2 .The panels correspond to different planet-to-star mass ratios as marked by the labels.Each planet is held in a fixed circular orbit with a p = 0.12 au.We overplot streamlines of the gas flow relative to the planet with white curves.The figure is also available as an online movie that scans all values of a p for which the static torque was measured.Corresponding torque measurements can be found in the top row of Fig. 12 (blue data points).TOI-216b), otherwise an extremely large α MRI would be needed to trigger the trap and that might be difficult to explain with our current understanding of viscosity-generating processes.For comparison purposes, let us recall that Flock et al. (2017) find α MRI 10 −1 and α DZ 10 −3 from MHD models of the inner disk turbulence.These values are consistent with the existence of the planet trap for TOI-216b and similar planets.
Future works can bring additional constraints in conjunction with our study.If, for example, future MHD models were to discover that viscosity transitions are rather wide than narrow, the bottom panel of our Fig. 8 would immediately imply α MRI 10 −1 for TOI-216b.
Finally, since the trapping becomes progressively difficult for q larger than that of TOI-216b, a question arises as to whether there could be a dynamical contribution to the Neptunian desert.We speculate that general properties of the viscosity transition at the inner rim might be such that sub-Neptunes are trapped easily, Neptunes like TOI-216b are trapped occasionally (as their abundance is quite small) and super-Neptunes are not trapped.Their continuing inward migration would thus contribute to their observational paucity -it would "clear the desert".Although the likelihood of this effect cannot be assessed at the moment, we merely wish to point it out as a possibility.If our speculations were true, additional questions would require answers, for instance, what happens with super-Neptunes once they reach the disk edge carved by the magneto-spheric cavity?Are these planets simply evaporated or could they stop once their outer Lindblad resonances migrate out of the disk?Could they become hot Jupiters by merging with companion super-Neptunes that have migrated in a similar fashion?

Relevance for N-body simulations with planet migration
Our work has important implications for planet population synthesis as well as for N-body simulations with prescribed migration.The latter cannot account for the effects shown in our study because the influence of the planet on the inner rim structure and the development of the front-rear asymmetry in the co-rotation region have not yet been included in the migration formulae (Paardekooper et al. 2011;Jiménez & Masset 2017;Kanagawa et al. 2018).Therefore, N-body simulations that lead to a formation of a sub-Neptune or Neptune as the innermost planet of the planet chain might suffer from systematic errors because such an innermost planet would feel an incorrect torque.Our study thus adds to the argument of Brasser et al. (2018) who demonstrate that N-body models with planet migration in the vicinity of the disk edge require refinements.
To highlight the incompleteness of torque formulae, we computed the normalized torque by applying the recipe of Paardekooper et al. (2011, see their Sect. 5.7) to the radial profiles of our non-isothermal disk (unperturbed by the planet).We chose α DZ = 10 −3 and α MRI = 5 × 10 −2 , which is the same disk as in Figs.11 and 13.Since the formulae of Paardekooper et al. (2011) are only applicable in the Type I regime of planetary migration, a criterion is needed to define gap-opening planets that migrate in the Type II regime and exclude them from the calculation.One such criterion is that of Crida et al. (2006), which remains to be used (e.g.Flock et al. 2019;Venturini et al. 2020b;Emsenhuber et al. 2021) even though updated criteria do exist (Kanagawa et al. 2018;Duffell 2020).
Combining Paardekooper et al. (2011) and Crida et al. (2006), we obtained the migration map shown in the top panel of Fig. 14.After comparing with the results of our simulations (horizontal lines in Fig. 14), it is clear that the map correctly predicts the trap location but it incorrectly overpredicts the migration trap efficiency -trapping predicted for planets as massive as M p 115 M ⊕ or q = 4.5 × 10 −4 (more massive than TOI-216b by a factor of six).This error is caused by the failure of the gapopening criterion of Crida et al. (2006) to account for shallow gaps and thus the formulae of Paardekooper et al. (2011) are applied even to planetary masses for which the disk response is already non-linear.
Next, we tested more recent migration formulae.Specifically, we calculated the Lindblad torque Γ L and the co-rotation torque A63, page 15 of 20 Fig. 14.Qualitative comparison between our hydrodynamic simulations (horizontal lines) and torque formulae (colour gradient) that are often used in N-body simulations of planet formation.Migration maps consisting of red and blue colour gradients mark the regions of the parametric space in which the torque formulae predict outward and inward migration, respectively.The black isocontour is where the predicted torque is exactly zero.Orange and cyan lines mark the extent of semimajor axes for which our simulations resulted in positive and negative torques, respectively (see Fig. 12).In the top panel, γΓ/Γ 0 is calculated following Paardekooper et al. (2011) for all planetary masses that do not violate the gap-opening criterion of Crida et al. (2006) (planets that do violate Crida's criterion belong to the grey region).In middle and bottom panels, the torque is calculated according to Jiménez & Masset (2017) and blended between Type I and Type II regimes using the recipe of Kanagawa et al. (2018).Top and middle panels are based on the unperturbed radial profiles of our non-isothermal disk with the viscosity transition α DZ = 10 −3 and α MRI = 5 × 10 −2 (see Fig. 11); viscosities in the bottom panel are α DZ = 2.5 × 10 −4 and α MRI = 5 × 10 −2 .We recall that the planet mass M p is obtained from q assuming the stellar mass of TOI-216, M = 0.77 M .
Γ C according to Jiménez & Masset (2017, see their Sect. 5.1) and then we blended them together as suggested by Kanagawa et al. (2018): where K is given by Eq. ( 2).The blending should ensure a smooth transition between Type I and Type II migration and the resulting torque should be applicable even when the planet opens a gap.
The result is shown in the middle panel of Fig. 14 and, at first glance, it seems to match our simulations very well -although the trap location is shifted inwards by 0.01 au compared to our simulations, the maximum planetary mass that can be trapped is found correctly.However, further experiments revealed that the match is rather coincidental.By repeating the calculation for the disk with α DZ = 2.5 × 10 −4 and α MRI = 5 × 10 −2 , as shown in the bottom panel of Fig. 14, we discover that the island of outward migration is much less sensitive to the viscosity contrast compared to our simulations.While we can barely trap the planet with q = 5 × 10 −5 in our simulations with a low-viscosity dead zone, the torque formulae predict a similar trap as in the middle panel of Fig. 14 and thus they overestimate the trapping efficiency yet again.Therefore, we conclude that the available torque formulae indeed fail to recover the results found in our hydrodynamic simulations.

Comparison to other works
Our results indicate a different regime of trapping compared to Masset et al. (2006), as we already pointed out in previous sections.In addition to that, let us emphasize that our mechanism differs from the vortex-aided trapping as well (Regály et al. 2013;Ataiee et al. 2014) because in that scenario, a Rossby vortex needs to be formed at a viscosity transition (or at a density bump) prior to the planet reaching the vortex location.The trap is then triggered during the orbital crossing of the planet and the vortex.In our case, instead, if a vortex forms, it does so gradually inside the co-rotation region of the planet.The vortex itself is not a necessary prerequisite for the trapping, it is rather a secondary result of gas concentration in the asymmetric co-rotation region.We also detect vortices forming outside the planet's co-rotation region but these are not long-lived and do not seem to contribute to trapping.
We also find differences compared to Romanova et al. ( 2019), likely because the authors did not account for viscosity transitions, their simulation timescales were short (did not allow for full gap opening), the whole planet mass was inserted immediately in the simulation and the planet was allowed to migrate straightaway.On the other hand, Romanova et al. ( 2019) used more realistic 3D simulations and it is left for future work to verify whether our findings can be recovered in 3D as well.
Unlike in Faure & Nelson (2016), our non-isothermal simulations do not form any RWI vortices at the transition between the dead and active zones that would undergo cyclic evolution and migration.We think that some of our model features would need to be modified in order to recover results of Faure & Nelson (2016), for example, the viscosity would need to be temperaturedependent (not radially dependent), our dead zone would need to have zero viscosity, or our aspect ratio would need to be larger.
A firm overlap can be found between our study and Hu et al. ( 2016) who explore the migration trap at the inner disk rim in the context of the inside-out planet formation of close-in A63, page 16 of 20 exoplanetary systems (Chatterjee & Tan 2014).By looking at the radial surface density profile perturbed by their most massive planet with q ≈ 3.5 × 10 −5 (see the top panel of their Fig.3), we see that it roughly corresponds to the perturbation produced by our least massive planet with q = 1.5 × 10 −5 (see panel a of our Fig.6).The slight difference in the gap-opening ability is likely caused by differences in the aspect ratio of the background disk model.However, one effect of Hu et al. (2016) that is not considered in our model is the outward shift of the dead-zone inner edge related to the penetration of ionizing X-rays into the gaped disk.
Finally, let us discuss the streamline topology.The appearance of a front island of librating streamlines is reminiscent of results of Ogilvie & Lubow (2006, see their Fig. 4) where the deformation of the co-rotation streamlines is a result of planet migration.The front-rear asymmetry is additionally reminiscent of Type III migration (Masset & Papaloizou 2003;Pepliński et al. 2008a), although that regime is usually considered in the context of runaway inward migration (cf.Pepliński et al. 2008b) and not in the context of planet traps.We point out that in our case, compared to studies mentioned above, the asymmetry of the co-rotation region is recovered regardless of planet migration, in both static and dynamical numerical experiments.Therefore, as stated in Paardekooper (2014), the formation of a one-sided mass deficit in the co-rotation region that we reported is likely a natural consequence of the planet finding itself on an initially steep density gradient (the planet is radially sandwiched between a low-density/high-viscosity inner region and a highdensity/low-viscosity outer region).For completeness, let us note that our streamline topology is azimuthally antisymmetric to the one found in McNally et al. (2017, e.g. their Figs. 4 and 6).

Limitations and future work
Due to the need to perform simulations over relatively long timescales, our disk model inevitably contains several simplifications and some of them should be tested in future to capture the physics of the inner disk rim in its entire complexity.The main two drawbacks are (i) the 2D approximation and (ii) the viscosity prescription.While it is relatively easy to convert our model to 3D (at least the locally isothermal one), improving the viscosity is not straightforward.One can at least think of a temperature-dependent viscosity, as already pointed out, but an equally important question is what happens with the disk turbulence in the presence of the planetary gap.If the gap becomes deep in the dead zone (before approaching the inner rim), can the gap centre become exposed and ionized to increase the turbulence level locally?And, since our viscous disk is laminar, what would be the evolution in a truly turbulent disk?Only future works can answer these questions.
Regarding the planet itself, we assumed that its mass remains constant.In other words, we neglected pebble accretion and gas accretion.While the omission of the former can be advocated because the pebble isolation mass is typically low in the inner disk (e.g.Bitsch 2019;Venturini et al. 2020b), gas accretion could still be ongoing.We decided not to include gas accretion in our global 2D models because of two obstacles -the implementation would have to rely on accretion rates adopted from other studies (e.g.Machida et al. 2010;Bodenheimer & Lissauer 2014) and the planet would have to act like a sink particle (while in an ideal situation, one would like to resolve a 3D gas flow all the way to the planetary envelope; e.g.Lambrechts & Lega 2017;Schulik et al. 2019).Nevertheless, gas accretion should be considered in future because it can change the gap profile (Bergez-Casalou et al. 2020) and thus affect the asymmetry of the co-rotation region as well.
As a side result, we find Type III migration for q = (1.5-4.5)× 10 −5 shortly before these planets reach the viscosity transition.The importance of this migration type for the assembly of close-in exoplanetary systems should thus be reassessed in future studies.

Conclusions
Motivated by recent claims that the innermost sub-Neptune-like exoplanets in multi-planetary systems cluster at 0.1 au (Mulders et al. 2018) and that such clustering is caused by the migration trap acting at a viscosity (or density) transition at the inner disk rim (Masset et al. 2006;Flock et al. 2019), we performed 2D locally isothermal and non-isothermal simulations to revisit the trap mechanism for sub-Neptunes and Neptunes (although our actual planet mass interval extended beyond these types of planets as well).Our working hypothesis was that close-in sub-Neptunes and Neptunes are actually able to open prominent gaps because of the very low aspect ratio h 0.02-0.03that is expected at the disk rim (Flock et al. 2019;Nesvorný et al. 2022).Consequently, the migration trap of Masset et al. (2006) might fail because it relies on the positive vortensity-driven corotation torque which might be suppressed when the co-rotation region becomes emptied due to a gap opening.As our reference representative of close-in Neptunes, we chose TOI-216b since its orbital properties and its firm 2:1 resonance with TOI-216c can be explained in a scenario that relies on the migration trap at the inner rim (Nesvorný et al. 2022).
We discover a new regime of the migration trap that is facilitated by the development of a front-rear asymmetry in the co-rotation (or horseshoe) region of the planet.The front horseshoe region develops an island of librating streamlines that accumulate a gas overdensity.Since the gas overdensity is leading the orbital motion of the planet, it exerts a positive torque and it can balance the negative torque of spiral wakes, thus halting the migration.In our locally isothermal runs, we find that the overdensity can turn into a Rossby vortex if the viscosity contrast is large and if q ≥ 6 × 10 −5 .The trap, however, persists even in the presence of the vortex.
For the parameters that we explored, the new regime of the migration trap operates roughly for planet-to-star mass ratios q (4-8) × 10 −5 .In this mass range, the planets open moderate gaps (for h 0.023 that we considered) and they locally erase the initial surface density transition, replacing it with a bump-like surface density distribution centred on the planet's corotation.For lower q, the gaps remain shallow and the trapping proceeds as suggested by Masset et al. (2006).For higher q, the total gas mass in the co-rotation region is further reduced and its net gravitational torque is no longer capable of sustaining the trap.
Instead of varying the planet mass, similar effects are found if the planet mass is fixed and the disk viscosity is varied instead.This is because the gap-opening ability increases with decreasing viscosity.By studying migration of TOI-216b across viscosity transitions with different dead-zone viscosities α DZ and MRI-active viscosities α MRI , we are able to find constraints for their combinations that promote the trapping.In summary, when α DZ = 10 −3 , TOI-216b can only be trapped for α MRI 5 × 10 −2 .Similarly, α DZ = 5 × 10 −4 requires α MRI 7.5 × 10 −2 for the trap to operate.Finally, the trap never occurs for α DZ ≤ 2.5 × 10 −4 , at least for the considered range of α MRI that had an upper A63, page 17 of 20 A&A 666, A63 (2022) limit of 10 −1 .Our results agree between the locally isothermal and non-isothermal disks, which makes them robust and independent of the disk thermodynamics, although the viability of some of our model assumptions (see Sect. 5) should be tested in future studies.Additionally, we find very similar locations of the migration trap when comparing dynamical simulations to static simulations with a fixed planetary orbit.In other words, the existence of the trap does not seem to strongly depend on the migration history of the planet.
Our findings imply that by studying sub-Neptunes and Neptunes whose formation requires a migration trap at the inner rim, it is possible to place constraints on the viscosity in central parts of their natal protoplanetary disks.We also speculate that there might be a dynamical contribution to the existence of the Neptunian desert because planets more massive than Neptunes are difficult to trap in general.Finally, we point out that N-body models with prescribed migration might suffer from systematic errors at the inner rim because the available torque formulae seem to overpredict the maximum planet mass that can be trapped (see Sect. 5.2).
To remove the dependency on H, we take H = c s √ γΩ = γP/ρ γGM /r 3 = (γ − 1) c V T i,out r 3 GM , (B.5) with the interior temperature T i,out = 2 −1/4 T s,out .The latter relation between the interior and surface temperature is based on the assumption of perfect splitting of the irradiating flux in halves: one half is radiated away and one half goes deeper into the disk (Dullemond et al. 2001).Finally, .

(B.6)
A generalized solution for the midplane temperature of an optically thick irradiated disk can be written as (see also Ueda et al. 2017) where the factor 1/2 again accounts for the flux splitting.

Appendix C: Supplementary figures
In this section, we present additional figures to support our claims related to Sect.3.4, Fig. 4 and the peculiar cases of migration reported therein.First, we pointed out that planets with q ≤ 4.5 × 10 −5 undergo runaway inward Type III migration before becoming trapped at the viscosity transition.Fig. C.1 shows the density distribution and streamlines near the planet with q = 1.5 × 10 −5 migrating through the dead zone with α DZ = 10 −3 and approaching the MRI-active zone with α MRI = 2.5 × 10 −2 .The displayed state corresponds to t mig = 300 P orb in bottom panel of Fig. 4 and we selected the case in which the Hill cut of the disk material for the torque evaluation was considered.We can see that there is a clear gas mass deficit in the front horseshoe region and thus the negative torque felt by the planet is enhanced.This explains why the planet migrates rapidly towards the disk edge.Fig. C.2 shows the case with q = 3.6 × 10 −4 (the other parameters remain the same) for which the planetary migration stagnated even before reaching the viscosity transition.The stagnation occurs due to the excited eccentricity of the planet, e 0.05 ∼ 2h, because the planet starts to encounter the gap walls (Duffell & Chiang 2015).However, this case should be approached with caution as it might suffer from numerical artefacts.The reason is that for the given combination of q, h and α DZ , the planet gets very efficient in perturbing the disk.Consequently, the outer part of the disk becomes eccentric first.The eccentric disk then starts to wobble around the domain origin (because the system barycentre becomes offset from the star) and, because we do account for the indirect terms of the gravitational potential, the wobble is sustained.But at the same time, the wobbling disk becomes inconsistent with our boundary conditions that are supplemented with wave-killing zones.These wave-killing zones essentially enforce axial symmetry close to the domain boundary that does not match the eccentric disk.We think that future work is needed in order to exclude the possibility of a nonphysical feedback loop between the eccentric disk, a gap-opening planet with a deep gap and our boundary conditions.

Fig. 1 .
Fig.1.Temporal evolution of the semi-major axis a p of migrating TOI-216b.Horizontal axis measures the migration time t mig (we do not show the simulation stages during which the planet was non-migrating).Solid lines are for α DZ = 10 −3 and dashed lines are for α DZ = 5 × 10 −4 .Different values of α MRI are distinguished by colour and specified in the plot legend.The shaded area marks the disk region over which α DZ increases to α MRI .Outwards from the shaded area, there is the dead zone.The MRI-active zone is located inwards from the shaded area.

FigFig. 2 .
Fig.2.Evolution of the radial surface density profile for four limit-case viscosity transitions.Simulations shown here combine α DZ = 10 −3 (left column) and 5 × 10 −4 (right column) with α MRI = 10 −1 (top row) and 10 −2 (bottom row).The black curve corresponds to the unperturbed disk.Coloured curves correspond to different simulation times t as indicated by the plot legend.The planet is kept fixed for t = 10 4 and 1.5 × 10 4 P orb when α DZ = 10 −3 and 5 × 10 −4 , respectively.At larger times, the planet is allowed to migrate freely.The vertical dashed line marks the final location of the planet in simulations where trapping occurs.We point out that by setting Ṁ based on α DZ (see Sect. 3.1), the unperturbed density profile in the dead zone is identical in each simulation.
Figure2shows the temporal evolution of the radial surface density profile Σ(r) for four selected simulations that combine the minimum and maximum α-viscosities from our parametric range.To obtain the Σ(r) profiles, we calculated the arithmetic mean from Σ(r, θ) in each radial ring of the domain.Figure 2 confirms our predictions from Nesvorný et al. (2022) and Sect. 1 -the gap in the dead zone indeed becomes quite deep, with the minimum values Σ/Σ 0 0.03 when α DZ = 10 −3 and Σ/Σ 0 0.008 when α DZ = 5 × 10 −4 .We also immediately see a substantial difference in the gap profile between simulations in which the planet is trapped and those in which it continues migrating inwards.For simulations where the trapping fails (bottom row in Fig.2), we see that the gap mostly retains its characteristic shape even when the planet is crossing the viscosity transition -the gap centre still appears to be close to an inner-outer symmetry, with a small central peak that reflects the mass accumulation in the Hill sphere of the planet.The gap walls do exhibit an inner-outer asymmetry, with the outer wall becoming steeper.As a result, the outer spiral arm contains more gas mass relative to the inner spiral arm and thus drives a strong negative Lindblad torque that is responsible for the episodes of fast inward migration in Fig.1.

Fig. 3 .
Fig. 3. Temporal evolution of the gas surface density Σ (shown in a logarithmic colour gradient) in a simulation with α DZ = 10 −3 and α MRI = 10 −1 .In each panel, the depicted portion of the disk is centred on TOI-216b.A rectangular projection of azimuthal coordinates is used to highlight disk features orbiting at the same radial separation.Streamlines of the gas flow are overlaid as white oriented curves.Individual panels are labelled with the migration time t mig (corresponding to Fig.1) and the instantaneous orbital radius of the planet r p .Panels a-d capture important phases of migration of TOI-216b towards the viscosity transition, panel e shows a state shortly after the migration stalls and panel f shows the final state of our simulation.The figure is also available as an online movie showing the temporal evolution between t mig = 16 000 and 17 000 P orb .

Figure 3
Figure3shows several snapshots from our hydrodynamic simulation of migrating TOI-216b that becomes trapped at the viscous transition between α DZ = 10 −3 and α MRI = 10 −1 .Panel a shows TOI-216b while it is still positioned outside the outer edge of the A63, page 7 of 20

Fig. 5 .
Fig.5.As in Fig.3, but now for different values of the planet-to-star mass ratio q as indicated by labels of individual panels.The viscosity transition is parameterized with α DZ = 10 −3 and α MRI = 5 × 10 −2 (see top panel of Fig.4for corresponding migration tracks of individual planets).In panels a-e that exhibit the migration trap, we show the state of the simulation at t mig = 40 000 P orb .In panel f , the planet is not trapped and thus we show the state at t mig = 9000 P orb when the planet is roughly at the same radial distance as in panel e.super-Earths or mini-Neptunes become trapped at the viscosity transition efficiently.We refer the reader to Appendix C and Fig. C.1 for an extended discussion.The most massive planet in the bottom panel of Fig.4(purple curve) stops migrating before reaching the viscosity transition because it develops a significant eccentricity equal to e = 0.05 and then the torque becomes modified by interactions with the gap walls(Duffell & Chiang 2015).However, it is possible that this simulation suffers from numerical artefacts and an extended study would be needed to exclude or confirm such a possibility.The reason for doubt is that the dead zone properties are identical for both panels of Fig.4yet the planet behaves differently in the bottom panel.We provide additional explanatory points in Appendix C and Fig. C.2.In Fig.5, we show 2D maps of the gas density and streamlines in the vicinity of each considered planet for α DZ = 10 −3 and α MRI = 5 × 10 −2 .With the increasing planetary mass, the gap-opening capability becomes more efficient and we see the following trends: (i) the surface density transition at the inner rim becomes gradually erased and replaced by a density bump centred in the co-rotation region; (ii) the density in the co-rotation region and thus the total mass of gas that can exert the co-rotation torque decreases on average (see also Fig.6); (iii) the streamline topology changes from regular horseshoes towards a dominant A63, page 9 of 20

Fig. 8 .
Fig. 8.As in Fig. 1, but for fixed α DZ = 10 −3 and different values of the viscosity transition width parameter ∆r MRI = H = 0.023 au (top) and ∆r MRI = 4H = 0.092 au (bottom).We point out that the total width of the viscosity transition (as marked by the shaded area) is not exactly equal to ∆r MRI ; the influence of the parameter is described by Eq. (10).

Fig. 9 .
Fig. 9. Locally isothermal simulations with a non-migrating planet.Top: Unperturbed surface density profile (solid lines, left vertical axis) across the viscosity transition (dashed lines, right vertical axis).Bottom: Measurements of the normalized static torque at various orbital radii (points) scanning the viscosity transition.In both panels, α DZ = 10 −3 and different values of α MRI are distinguished by colours and symbols (see the plot legend).

Fig. 10 .
Fig. 10.Surface density and streamlines for viscosities α DZ = 10 −3 , α MRI = 10 −1 and TOI-216b in a fixed circular orbit at a p = 1.1 au.The plot serves as a comparison with panel f of Fig. 3 where the planet was migrating.The figure is also available as an online movie showing the variations in Σ and streamlines for all values of a p considered in our static torque measurements.

Fig. 11 .
Fig. 11.Radial profiles of the surface density Σ (top, left vertical axis, solid curves), α viscosity (top, right vertical axis, dashed curves), temperature T (bottom, left vertical axis, solid curves) and aspect ratio h (bottom, right vertical axis, dashed curves).A comparison is shown between our 2D non-isothermal model (black curves) and a 3D vertically resolved radiative hydrostatic model (red curves).The bottom panel is split into three regions discussed in Sect.2.2, namely the dust halo (B), the condensation front (C) and the flared disk (D) (e.g.Ueda et al. 2017).
0.110 0.115 0.120 0.125 0.130 semi-major axis a p [au] 0.110 0.115 0.120 0.125 0.130 semi-major axis a p [au] Fig. C.1.Surface density and streamlines for a migrating super-Earthsized planet with q = 1.5 × 10 −5 at t mig = 300 P orb after the planet was released.The corresponding migration track is shown in bottom panel of Fig. 4 (α DZ = 10 −3 and α MRI = 2.5 × 10 −2 ).There is an excess mass in the rear horseshoe region and the planet thus migrates in the runaway Type III regime.

Fig
Fig. C.2. Surface density for a migrating planet with q = 3.6 × 10 −4 shortly before the end of the simulation.The corresponding migration track is shown in the bottom panel of Fig. 4. The planet is strongly eccentric (e 0.05), which manifests itself via the eccentricity of the gap and the relative position of the planet with respect to gap edges.We point out that this unusual result (compared to our remaining simulations) might be caused by a numerical artefact.

Table 1 .
Parameters varied in our simulations.

Table 2 .
Fixed parameters of the locally isothermal simulations.

Table 3 .
Fixed parameters of the non-isothermal simulations.