© O. Chrenko et al. 2022

1 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 4R_⊕; 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 often 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¹ (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 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, 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 Neptune-class TOI-216b (M_p = 0.059 M_Jup, P_orb ≃ 17.1 d), and an outer half-Jupiter TOI-216c (Mp = 0.56 M_Jup, Porb ≃ 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 (1)

where Σ_gap is the minimal surface density at the gap bottom, Σ₀ is the unperturbed surface density, and (2)

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 disk². Before reaching the inner rim, the planet is embedded in the dead zone where α = 10⁻⁴−10⁻³ 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⁻⁵ (corresponding to a ten-Earth-mass planet orbiting a Solar-mass star), one obtains³ Σ_gap/Σ₀ ≃ 0.018−0.15. Such an excavation of the co-rotation 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⁻⁵ or 4.5 × 10⁻⁵, 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⁻⁵, 10⁻⁴) 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 & Konigl 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.

2 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 protoplan-etary 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.

2.1 Locally Isothermal 2D Model

Governing equations of the model read as follows: (3) (4)

where Σ is the gas surface density, υ is the velocity vector, P is the pressure, ⊤ is the viscous stress tensor, and ϕ is the gravitational potential. The equation of state is given by (5)

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 (6)

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 under-resolved 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): (7)

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 (8)

where the kinematic viscosity v was calculated using the alpha prescription (Shakura & Sunyaev 1973) (9)

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) (10)

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 (11)

and the azimuthal velocities υ_θ were set by the centrifugal balance between the gravity of the star and the pressure support within the disk. The initial conditions for Σ, υ_r and υ_θ were also used as boundary conditions. To provide a smooth transition between the boundaries and the perturbed disk, we damped Σ and υ_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).

2.2 Non-isothermal 2D Model

We relaxed the locally isothermal approximation from previous Sect. 2.1 by replacing Eq. (5) with (12)

where γ is the adiabatic index and ℰ is the internal energy density. To describe the temporal evolution of ℰ, 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) (13)

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∇ · υ accounts for the compres-sional 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 (14)

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 ℰ = ρc_VT, 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) and D'Angelo & Marzari (2012): (15)

with (16)

In writing Eq. (16), we introduced a single grey opacity κ, which was considered uniform.

Ultimately, we aimed for our model to be applicable at the inner rim of protoplanetary disks, which is thermodynamically complex. Flock et al. (2016, 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 D⁴ because they contain the first possible trap⁵ 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 (17)

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 (18)

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 (19)

where (20)

and (21)

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: (22)

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 (23)

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 ℰ = ρc_VT. 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: (24)

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 v 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⁻⁵. 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 hydrody- namic 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 cells⁶. Next, we used T in the ghost cells to find v. Afterwards, υ_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 = Σ_av_a/v_g, where the subscripts 'g' and 'a' stand for the ghost and active cells, respectively. At the outer boundary, we used ∑_g = Ṁ/(3πv_g) in order to replenish the disk material in a uniform way. Finally, ℰ 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 υ_r so that Σ could adjust and reach an equilibrium. The target value of υ_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 υ_r towards the equilibrium disk state found at the end of the relaxation stage.

3 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⁻⁵. 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.

Fig. 1 Temporal evolution of the semi-major axis ap of migrating TOI- 216b. Horizontal axis measures the migration time tmig (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.

3.1 The Role of the a-viscosity Contrast for Trapping TOI-216b

We start by placing TOI-216b at a_p(t₀) = 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⁻⁴ and 10⁻³. For each of these α_DZ values, we vary the MRI-active viscosity as α_DZ = (10⁻², 2.5 × 10⁻², 5 × 10⁻², 7.5 × 10⁻², 10⁻¹). 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⁻⁸ M_⊙ yr⁻¹ when α_DZ = 10⁻³ and Ṁ = 5 × 10⁻⁹ M₀ yr⁻¹ when α_DZ = 5 × 10⁻⁴. 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⁻³ and 5 × 10⁻⁴, respectively. Only after a gap was formed and its opening slowed down did we let the planet migrate for 50000 P_orb.

Since a_p(t₀) = 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 (Durmann & 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 ap after the release of TOI-261b is shown in Fig. 1 for all studied viscosity 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.3 au, 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⁻³, the planet is trapped at the viscosity transition for α_MRI ≥ 5 × 10⁻², otherwise it continues migrating inwards. When α_DZ = 5 × 10⁻⁴, the trapping is only possible for α_MRI ≥ 7.5 × 10⁻². 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.

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 = 104 and 1.5 × 104 Porb 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.

3.2 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. 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.03 when α_DZ = 10⁻³ and Σ/Σ₀ ≃ 0.008 when α_DZ = 5 × 10⁻⁴.

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. 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 tmig (corresponding to Fig. 1) and the instantaneous orbital radius of the planet rp. 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 tmig = 16 000 and 17 000 Porb.

3.3 Surface Density and Gas Flow Evolution in the co-rotation Region

Figure 3 shows several snapshots from our hydrodynamic simulation of migrating TOI-216b that becomes trapped at the viscous transition between α_DZ = 10⁻³ and α_MRI = 10⁻¹. Panel a shows TOI-216b while it is still positioned outside the outer edge of the 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 inward- migrating 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 overden- sity 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.

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.

3.4 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⁷ q = (1.5, 3, 4.5, 6, 9, 18, 36) × 10⁻⁵. Fig. 4 shows the temporal evolution of the planetary semi-major axes a_p(t) for fixed α_DZ = 10⁻³ and two values of the MRI-active viscosity α_MRI = 2.5 × 10⁻² and α_MRI = 5 × 10⁻². 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 40000 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⁻⁵ 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 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. 4 yet 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⁻³ and α_MRI = 5 × 10⁻². 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 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⁻⁵.

In Fig. 5, planets with q ≤ 7.3 × 10⁻⁵ 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⁻⁵ 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⁻⁵−7.3 × 10⁻⁵) and finally to untrappable planets that clear their co-rotation regions too efficiently (e.g. q = 9 × 10⁻⁵).

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⁻⁵ and for these planets, it produces the positive torque necessary for their trapping. For q = 9 × 10⁻⁵, 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 (25)

where the vorticity component is perpendicular to the disk plane and the term related to the orbital frequency translates ζ⁻¹ to the inertial frame. The inverse potential vorticity has several diagnostic applications in our case. First, we look for the local maxima of ζ⁻¹ to trace vortices that are the usual outcomes of non-linear RWI growth⁸. Second, the vortensity-driven horseshoe drag scales as Γ_HS ∝ Σd log ζ⁻¹ /d log r (Masset et al. 2006) so that steep radial gradients of ζ⁻¹ 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⁻⁵. For the lowest planet mass, we see that a steep positive gradient of ζ⁻¹ 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 ζ⁻¹ 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 ∇ × υ and mixing with a lower Σ) reduce ζ⁻¹ along the outer separatrix of the rear horseshoe region.

The middle and bottom panels of Fig. 7 show that the gradient of ζ⁻¹ across the co-rotation is reduced with increasing planet mass, as we already anticipated. Moreover, we notice that a steep isolated local maximum of ζ⁻¹ only appears for q = 7.3 × 10⁻⁵. For q = 4.5 × 10⁻⁵, the peak values do not exceed those found for q = 1.5 × 10⁻⁵ and thus the planet with q = 4.5 × 10⁻⁵ 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⁻⁵ 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⁻⁵ than for q = 7.3 × 10⁻⁵ (panels c and e, respectively, 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⁻⁵ 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).

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. 4 for corresponding migration tracks of individual planets). In panels a-e that exhibit the migration trap, we show the state of the simulation at tmig = 40000 Porb. In panel f, the planet is not trapped and thus we show the state at tmig = 9000 Porb when the planet is roughly at the same radial distance as in panel e.

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 torb = 40000 Porb for trapped planets and at torb = 9000 Porb 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.

Fig. 7 2D distribution of the inverse potential vorticity ζ−1 for the same simulation set as in Figs. 5 and 6 and for q = 1.5 × 10−5 (top), 4.5 × 10−5 (middle) and 7.3 × 10−5 (bottom). Vertical dashed lines mark the approximate extent of the horseshoe region following Jiménez & Masset (2017, their Eq. (30)). The colour scale is saturated to highlight features in the co-rotation region.

Fig. 8 As in Fig. 1, but for fixed αDZ = 10−3 and different values of the viscosity transition width parameter ∆rMRI = H = 0.023 au (top) and ∆rMRI = 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 ∆rMRI; the influence of the parameter is described by Eq. (10).

3.5 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⁻³, α_MRI = (10⁻², 2.5 × 10⁻², 5 × 10⁻², 7.5 × 10⁻², 10⁻¹) 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 30000 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⁻² compared to Fig. 1. If, on the other hand, the viscosity transition becomes flattened, the trapping becomes increasingly difficult and only occurs for the most turbulent MRI zone with α_MRI = 10⁻².

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 magne- tohydrodynamic studies should therefore improve the ionization profile and disk chemistry to carefully reassess the radial disk range over which the MRI is triggered.

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).

3.6 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 front- rear 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⁻³ and α_MRI = (2.5 × 10⁻², 5 × 10⁻², 10⁻¹). Additionally, the planetary semi-major axis remains fixed and we take a_p distributed evenly between 0.975 au and 1.15au 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 20000 P_orb to allow the disk torque to converge.

Figure 9 shows our measurements of the static torque Γ normalized to (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 in Fig. 9, outward migration is expected and vice versa. For α_MRI = 10⁻¹, 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⁻², 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×l0⁻²,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⁻³, α_MRI = 10⁻¹ 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.

Fig. 10 Surface density and streamlines for viscosities αDZ = 10−3, αMRI = 10−1 and TOI-216b in a fixed circular orbit at ap = 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 ap considered in our static torque measurements.

4 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 remaining fixed parameters of our non-isothermal simulations are summarized in Table 3.

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).

4.1 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 . The reference 3D simulation in Fig. 11 is based on the radiative hydrostatic approach of Flock et al. (2016, 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 (26)

To model the radius-dependent α-transition in the 2D simulation (see Eq. 10), we chose r_MRI = 0.11au 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.

4.2 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⁻⁵), 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⁻⁴, four values of α_MRI = (2.5, 5, 7.5, 10) × 10⁻² and three values of q = (5, 7.3, 9) × 10⁻⁵. 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 v 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 γΓ/Γ₀. 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 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⁻¹.

Looking at the case of q = 5 × 10⁻⁵, 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⁻⁴ and α_MRI = 2.5 × 10⁻². 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⁻⁴.

Focussing now on TOI-216b with q = 7.3 × 10⁻⁵, we see that α_DZ = 10⁻³ requires α_MRI ≲ 5 × 10⁻² in the MRI-active zone to sustain the trap and α_DZ = 5 × 10⁻⁴ leads to trapping only if α_MRI ≳ 7.5 × 10⁻². For α_DZ = 2.5 × 10⁻⁴, 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.

For the most massive planet shown in Fig. 12 with q = 9 × 10⁻⁵, only α_DZ = 10⁻³ and α_MRI ≳ 7.5 × 10⁻² trap the planet at the inner rim. For α_DZ ≤ 5 × 10⁻⁴, the trapping is not possible within our parametric range.

Finally, we investigate the surface density distribution and streamlines. For a comparison with non-isothermal simulations, we select the viscosity combination of α_DZ = 10⁻³ and α_MRI = 5 × 10⁻². 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⁻⁵, small and positive for TOI-216b, small and negative for q = 9 × 10⁻⁵).

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 region - a 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 overden- sity located in the front horseshoe region, in accordance with our locally isothermal simulations.

Fig. 12 Normalized static torque γΓ/Γ0 measured in our non-isothermal simulations as a function of the planetary orbital distance ap. 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 aDZ 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.

5 Discussion

5.1 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⁻² when α_DZ = 10⁻³ or α_MRI ≳ 7.5 × 10⁻² when α_DZ = 5 × 10⁻⁴) 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⁻⁴ for 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⁻¹ and α_DZ ≃ 10⁻³ 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⁻¹ 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?

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 ap = 0.12au. 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 ap for which the static torque was measured. Corresponding torque measurements can be found in the top row of Fig. 12 (blue data points).

5.2 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⁻³ and α_MRI = 5 × 10⁻², 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; Emsenhu- ber 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 is predicted for planets as massive as M_p ≃ 115 M_⊕ or q = 4.5 × 10⁻⁴ (more massive than TOI-216b by a factor of six). This error is caused by the failure of the gap- opening 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 Γ_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): (27)

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⁻⁴ and α_MRI = 5 × 10⁻², 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⁻⁵ 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.

Fig. 14 Qualitative comparison between our hydrodynamic simulations (horizontal lines) and torque formulae (colour gradient) that are often used in TV-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 semi- major 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 & Mas- set (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 Mp is obtained from q assuming the stellar mass of TOI-216, M* = 0.77 M⊙.

5.3 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 temperature- dependent (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 exoplanetary systems (Chatterjee & Tan 2014). By looking at the radial surface density profile perturbed by their most massive planet with q ≈ 3.5 × 10⁻⁵ (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⁻⁵ (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 high- density/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).

5.4 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⁻⁵ 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.

6 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.03 that 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 co- rotation 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⁻⁵. 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. 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 co-rotation. 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⁻³, TOI-216b can only be trapped for α_MRI ≳ 5 × 10⁻². Similarly, α_DZ = 5 × 10⁻⁴ requires α_MRI ≳ 7.5 × 10⁻² for the trap to operate. Finally, the trap never occurs for α_DZ ≤ 2.5 × 10⁻⁴, at least for the considered range of α_MRI that had an upper limit of 10⁻¹. 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 Nep- tunes 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).

Acknowledgements

This work was supported by the Czech Science Foundation (grant 21-23067M) and the Ministry of Education, Youth and Sports of the Czech Republic through the e-INFRA CZ (ID:90140, LM2018140). The work of O.C. was supported by the Charles University Research program (No. UNCE/SCI/023). M.F. was supported by the European Research Council (ERC) project under the European Union's Horizon 2020 research and innovation programme number 757957. We wish to thank an anonymous referee whose valuable comments allowed us to improve this paper.

Appendix A Viscous Heating Term

The viscous dissipation term in full 3D spherical coordinates is (e.g. Mihalas & Weibel Mihalas 1984; D'Angelo & Bodenheimer 2013) (A.1)

where τ_ij are the components of the viscous stress tensor. Using vertical averaging and omitting zero terms in 2D, we obtain (A.2)

The remaining vertical component is non-zero because in (A.3)

only the strain tensor component D_ϕϕ is zero. Thus the final form for the viscous heating term is (A.4)

Appendix B Stellar Irradiation Term

Irradiating rays propagate from the star and they are impinging the τ ≈ 1 surface of the disk with the total flux of (B.1)

where α is the grazing angle and L_* is the stellar luminosity. Eq. (B.1) is applicable when the whole star is visible from an arbitrary point of the irradiated surface. This is consistent with the assumption of an optically thin region inwards from the condensation front of dust grains⁹ (Dullemond et al. 2001).

We assume that the τ ≈ 1 surface is caused by dust grains and we write for the energy balance of a dust grain (B.2)

where H_ph is the photospheric height above the midplane, C_bw = 4 is the backwarming factor (Kama et al. 2009), and we simplified the grazing angle prescription (cf. Chiang et al. 2001). Close to the star, the first geometrical term in the square brackets dominates and, using , we obtain the surface temperature (B.3)

Farther out, the second geometrical term in the square brackets prevails. To find a solution for the surface temperature, we assume that the disk obeys the flaring disk principle with the equilibrium slope d lnH_ph/d ln r = 9/7 (Chiang & Goldreich 1997). The photospheric height is by a factor of few larger than the pressure scale height, H_ph = f_phH, and this leads to (B.4)

To remove the dependency on H, we take (B.5)

with the interior temperature T_i,out = 2^−1/4T_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) (B.7)

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⁻⁵ 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⁻⁵ migrating through the dead zone with α_DZ = 10⁻³ and approaching the MRI-active zone with α_MRI = 2.5 × 10⁻². 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⁻⁴ (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. C.1 Surface density and streamlines for a migrating super-Earth- sized planet with q = 1.5 × 10−5 at tmig = 300Porb 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. 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.

References

All Tables

All Figures

Fig. 1 Temporal evolution of the semi-major axis ap of migrating TOI- 216b. Horizontal axis measures the migration time tmig (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.

In the text

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 = 104 and 1.5 × 104 Porb 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.

In the text

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 tmig (corresponding to Fig. 1) and the instantaneous orbital radius of the planet rp. 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 tmig = 16 000 and 17 000 Porb.

In the text

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.

In the text

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. 4 for corresponding migration tracks of individual planets). In panels a-e that exhibit the migration trap, we show the state of the simulation at tmig = 40000 Porb. In panel f, the planet is not trapped and thus we show the state at tmig = 9000 Porb when the planet is roughly at the same radial distance as in panel e.

In the text

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 torb = 40000 Porb for trapped planets and at torb = 9000 Porb 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.

In the text

	Fig. 7 2D distribution of the inverse potential vorticity ζ−1 for the same simulation set as in Figs. 5 and 6 and for q = 1.5 × 10−5 (top), 4.5 × 10−5 (middle) and 7.3 × 10−5 (bottom). Vertical dashed lines mark the approximate extent of the horseshoe region following Jiménez & Masset (2017, their Eq. (30)). The colour scale is saturated to highlight features in the co-rotation region.
In the text

	Fig. 8 As in Fig. 1, but for fixed αDZ = 10−3 and different values of the viscosity transition width parameter ∆rMRI = H = 0.023 au (top) and ∆rMRI = 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 ∆rMRI; the influence of the parameter is described by Eq. (10).
In the text

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).

In the text

	Fig. 10 Surface density and streamlines for viscosities αDZ = 10−3, αMRI = 10−1 and TOI-216b in a fixed circular orbit at ap = 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 ap considered in our static torque measurements.
In the text

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).

In the text

Fig. 12 Normalized static torque γΓ/Γ0 measured in our non-isothermal simulations as a function of the planetary orbital distance ap. 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 aDZ 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.

In the text

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 ap = 0.12au. 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 ap for which the static torque was measured. Corresponding torque measurements can be found in the top row of Fig. 12 (blue data points).

In the text

Fig. 14 Qualitative comparison between our hydrodynamic simulations (horizontal lines) and torque formulae (colour gradient) that are often used in TV-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 semi- major 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 & Mas- set (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 Mp is obtained from q assuming the stellar mass of TOI-216, M* = 0.77 M⊙.

In the text

	Fig. C.1 Surface density and streamlines for a migrating super-Earth- sized planet with q = 1.5 × 10−5 at tmig = 300Porb 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.
In the text

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.

In the text