© ESO 2020

1 Introduction

Giant planets form in protoplanetary disks and their early evolution is driven by gravitational planet-disk interactions. Once a forming giant planet exceeds a certain critical mass, it starts to deposit enough angular momentum in the surrounding disk to overcome the viscous spreading of gas; the gas is expelled away from the planetary orbit and a gap is opened (Crida et al. 2006; Kanagawa et al. 2015). In response, the angular momentum exchange with the disk forces the planet to migrate. The migration regime in the presence of the depleted corotation region of the planet is referred to as Type II migration (Lin & Papaloizou 1986a,b).

The classical paradigm of Type II migration states that after the gap opening, the gas flow across the gap is blocked and the planet behaves as if it is frozen in the gap centre (Lin & Papaloizou 1986a,b). Since a large portion of the disk accretes onto the central star, the planet is radially displaced along with the accretion flow. According to the theory of the viscous evolution of razor-thin accretion disks (Shakura & Sunyaev 1973; Lynden-Bell & Pringle 1974; Frank et al. 2002), the accretion velocity and therefore the expected Type II migration velocity is (1)

where ν is the kinematic viscosity of the disk and r is the radial distance.

Recently, the classical paradigm has been challenged. Lubow & D’Angelo (2006) and Duffell et al. (2014) found that there are gas flows across the gap and the assumption of the classical paradigm is therefore rendered invalid. On top of that, Hasegawa & Ida (2013) argued that there is no valid physical principle that would fix the planet in the centre of the gap. Dürmann & Kley (2015)demonstrated that the Type II migration rate scales with the planet mass and the local disk mass rather than with the disk accretion flow.

Three novel views of the Type II migration mechanism have been formulated on the basis of 2D locally isothermal simulations. First, Kanagawa et al. (2018) argued that the Type II migration rate can be actually predicted using the Type I migration physics while taking into account the decreased gas density in the gap. Second, Robert et al. (2018) suggested that after the gap opening, the planet starts to migrate as dictated by the differential Lindblad torque of the spiral arms, which is usually negative (Ward 1986). As the planet migrates inwards, the gap has to follow but the diskreshapes with a certain lag, which is determined by the timescale of viscous spreading. Through this lag (since the disk has to adjust to displacements of the planet), Type II migration remains dependent on ν but the migration speed tends to be faster than v_r,visc unless the disk mass becomes small enough (Dürmann & Kley 2015). Third, Scardoni et al. (2020) performed long-term simulations and found that the planet indeed initially migrates faster than v_r,visc (in accordancewith Dürmann & Kley 2015; Robert et al. 2018). However, if the migration is allowed to proceed for about ~ 10³ orbital timescales, the drift rate eventually converges to v_r,visc. Nevertheless, the portion of the disk that the planet crosses before the migration slows down is substantial – the planet starting at r_p ends up at ≃ 0.2r_p.

Whatever the true physical mechanism, the aforementioned studies generally agree that if giant planets formed at several au (as required by the core accretion scenario; Pollack et al. 1996), their inward Type II migration would operate on a timescale shorter than the typical lifetime of protoplanetary disks (Nelson et al. 2000; Hasegawa & Ida 2013). In other words, a great number of giant planets would be lost and the survivors would likely become hot Jupiters. Such an outcome, however, would be inconsistent with observations, which have revealed that the majority of giant planets orbit at separations ≳ 1 au (Mayor et al. 2011; Cassan et al. 2012; Fressin et al. 2013; Santerne et al. 2016). The inconsistency between the theory and observations can only be alleviated if giant planets form at fairly large separations of ≃15–30 au by efficient accretion processes (Coleman & Nelson 2014, 2016; Bitsch et al. 2019; Johansen et al. 2019) or if there is a special mechanism that can slow Type II migration down (e.g. Kanagawa 2019).

An intriguing mechanism that we highlight here was suggested by Hallam & Paardekooper (2018). They proposed that Type II migration can be slowed down or even reversed when stellar irradiation by the central star is taken into account. The slowdown should work as follows. After the gap opening, the outer edge of the gap receives an increased amount of stellar irradiation, it becomes hotter, and puffs up. The vertical expansion of the disk boosts the local aspect ratio h = H∕r where the pressure scale height is (2)

where c_s is the sound speed, γ is the adiabatic index, and Ω is the local Keplerian frequency. Since the one-sided Lindblad torque scales as Γ_OS ~ Ch⁻³ (Ward 1997; Papaloizou et al. 2007), where C is a positive constant for the inner disk and a negative constant for the outer one, the increase of h at the heated outer edge reduces the outer negative Lindblad torque. The total Type II torque thus becomes reduced and can even become positive if the heating of the outer gap edge is sufficiently strong.

Although promising, the results of Hallam & Paardekooper (2018) were obtained in a 2D model with simplified thermodynamics. They only accounted for compressional heating and local cooling parametrized by a thermal relaxation timescale. The stellar heating of the outer gap edge was not directly modelled; its influence was mimicked by an ad hoc Gaussian perturbation of the radial profile of c_s.

The central aim of our paper is to explore the mechanism proposed by Hallam & Paardekooper (2018) in a full 3D radiation hydrodynamics (RHD) model with stellar irradiation. We focus on passively heated disks (Chiang & Goldreich 1997) with a constant accretion rate provided by the α-viscosity approximation (Shakura & Sunyaev 1973). We consider a single embedded gap-opening planet. We explore the heating of the outer gap edge and we perform measurements of the static torque in which the planet is held at a fixed radial distance and the torque is computed from the gas density field.

Additionally, we stress that performing a new 3D radiative study of Type II migration is worthwhile. The majority of recent advances have been acquired through 2D locally isothermal simulations (e.g. Dürmann & Kley 2015; Robert et al. 2018; Scardoni et al. 2020), simply because 3D simulations are numerically demanding and also because Kley et al. (2001), Bitsch & Kley (2010), and Fung & Chiang (2016) identified only relatively small differences between 2D and 3D. Our results show that the inclusion of stellar irradiation, which is inherently a 3D phenomenon, can modify Type II migration.

2 Model

2.1 Physical principles

The gas disk is modelled as a viscous non-isothermal continuum on an annulus in spherical coordinates (comprising of the radius r, azimuth θ, and colatitude ϕ). The disk gravitationally interacts with two point mass objects M_⋆ and M_p, which represent a central star and a single embedded planet, respectively.

Our numerical experiments are conducted using the hydrodynamic code FARGO3D (Benítez-Llambay & Masset 2016) extended with our implementation of the radiation physics (Chrenko & Lambrechts 2019). The set of partial differential equations that describe the evolution of the gas disk and the radiation field reads

where ρ denotes the volume density, t the time, v the flow velocity vector, P the pressure, Φ the gravitational potential of the star and the planet, T the viscous stress tensor, r the radius vector, ɛ the internal energy of the gas, κ_P the Planck opacity, σ the Stefan–Boltzmann constant, T the gas temperature, c the speed of light, Q_irr the irradiation heating term, E_R the radiative energy, and F the radiation flux. The state equation of ideal gas together with the flux-limited diffusion approximation (FLD; Levermore & Pomraning 1981; Kley 1989) are used as closure relations for the system of Eqs. (3)–(6) (see Chrenko & Lambrechts 2019).

Our aim is to model the thermal balance of the disk in the regime of passive heating. For this reason, the gas energy Eq. (5) only contains the compressional heating term − P∇⋅v and the stellar-irradiation term Q_irr while the heating by viscous dissipation is not considered¹. Although this might seem unrealistic, such a setup allows us to isolate the influence of stellar irradiation on Type II migration more easily. If viscous heating were considered, it would induce bumps in the radial profile of the aspect ratio h(r) (Bitsch et al. 2013) and these could cause undesirable self-shadowing effects because the optical surface of the disk (with respect to stellar photons) would become bumpy as well. Keeping Q_irr only, the disk assumes a flared profile with h ∝ r^2∕7 (Chiang & Goldreich 1997) and self-shadowing can only occur once the planet is inserted and starts to perturb the gas distribution.

Stellar irradiation is implemented following Dobbs-Dixon et al. (2010), Bitsch et al. (2013), and Kolb et al. (2013). The central star with a physical radius R_⋆ and an effective temperature T_⋆ represents a point radiation source with a luminosity of . It is assumed that stellar photons impinging on the disk propagate along radial rays only, following paths of constant azimuth and colatitude. The optical depth to the irradiating flux is integrated along the radial rays as (7)

where τ₀ is the optical depth at the inner radial boundary of the domain r_min and κ_⋆ is the disk opacity to stellar photons (the Planck opacity at T_⋆). The local heating due to stellar irradiation is given simply by an exponential attenuation of the flux as (8)

where dτ_⋆ is the increment of the optical depth acrossa grid cell of interest, S_cell is the irradiated cross-section of the cell, and V_cell is its volume.

Regarding the disk opacity, we assume it is dominated by sub-micron dust grains that trace the distribution of the gas density. We adopt a simple opacity model from Flock et al. (2019) and set κ_P =κ_R ≡ κ_dust = 700 cm² g⁻¹ where κ_R is the Rosseland mean opacity (which enters the model through the FLD). The opacity to stellar irradiation is κ_⋆ = 1300 cm² g⁻¹. To evaluate theopacity of a gas–dust mixture, each opacity value is scaled by the dust-to-gas ratio, which we choose as Z = 0.001. The somewhat lower value of Z reflects the fact that sub-micron grains tend to get depleted as the mass spectrum of the coagulation-fragmentation equilibrium peaks at large grains (Birnstiel et al. 2012; Flock et al. 2019). We check for the evaporation of dust grains by calculating the evaporation temperature according to Isella & Natta (2005). Dust-free regions would have a reduced opacity κ_gas ≃ 10⁻⁵ cm² g⁻¹ (Flock et al. 2019) but since we focus on disk regions further out from the inner disk rim, the local temperature is always below T_evap. Our simple opacity treatment is again motivated by our effort to keep the flared disk profile monotonic; any temperature-dependent opacity transition would change the local cooling rate and the aspect ratio profile would become more complex (e.g. Bitsch et al. 2013).

To mimic the disk accretion due to the angular momentum transport, we use the classical parametrization by the α viscosity (Shakura & Sunyaev 1973) (9)

where γ is the adiabatic index (the ratio of specific heats) and c_V is the specific heat at constant volume, and we used the ideal gas state equation to expand the right-hand side.

The gravitational potential generated by the star and the planet is (10)

where d is the cell-planet distance smoothed by the cubic spline f_sm of Klahr & Kley (2006). We use the characteristic smoothing length r_sm = 0.5 R_H where R_H is the Hill sphere radius of the planet. In our calculations, the self-gravity of the gas is neglected. Therefore, our model cannot correctly account for interactions between the circumplanetary and protoplanetary disks and as a correction, we exclude the inner region of the Hill sphere when evaluating the disk-driven torque. The cut-off function is (Crida et al. 2008) (11)

and we use p = 0.6 (Robert et al. 2018).

The two-body star-planet interaction is computed with the IAS15 integrator (Rein & Spiegel 2015) from the REBOUND package (Rein & Liu 2012) and all simulations are performed in a reference frame corotating with the planet. The units used in the code are such that M_⊙ is the unit mass and au is the unit length. Furthermore, the gravitational constant as well as the ideal gas constant divided by the mean molecular weight are also equal to unity: G = 1; .

2.2 Parameters

The parameters that we use in our fiducial simulation are listed in Table 1. We focus on solar-type protostars (Baraffe et al. 1998; White et al. 2007) and our choice of R_⋆ and T_⋆ implies L_⋆ ≃ 1 L_⊙. Our fiducial planet is an analogue of a fully formed Jupiter. The grid is designed for consistency with previous 2D studies (e.g. Dürmann & Kley 2015; Robert et al. 2018) and resolves the R_H of a Jupiter-mass planet with approximately seven cells in each dimension.

Regarding the disk itself, our aim is to study accreting disks and we use the usual parametrization by the radial mass flux Ṁ together with α. The fiducial value of Ṁ = 10⁻⁸ M_⊙ yr⁻¹ is close to the mean value observed around solar-mass protostars in Lupus and Chamaeleon I (Manara et al. 2017; Mulders et al. 2017). The fiducial viscosity α = 10⁻³ is applicableto disks in which the angular momentum transport is facilitated by the hydrodynamic turbulence or disk winds (α ≃ 10⁻³–10⁻⁴; Nelson et al. 2013; Klahr & Hubbard 2014; Béthune et al. 2017) rather than magnetorotational instability (α ≃ 10⁻¹–10⁻²; Fromang & Nelson 2006; Flock et al. 2017). The latter is typically inactive at the considered location of the planet, a_p = 5.2 au (e.g. Matsumura & Pudritz 2005; Terquem 2008; Dzyurkevich et al. 2013). Viscosity transitions are omitted in our study for simplicity.

2.3 Simulation stages

Before simulating planet-disk interactions, it is necessary to ensure that the disk is in a thermal equilibrium and that its structure remains stationary over many dynamical timescales. Finding such an equilibrium state is not a straightforward task for an accreting non-isothermal disk because there is no generally valid analytic prescription for the disk profile. The difficulty lies in the following interplay. The disk mass at a given radius depends on Ṁ and ν. But ν itself is a function of T (Eq. (9)) and thus it is directly affected by the radiation reprocessing. The reprocessing, however, connects back to the gas distribution, which sets the relevant optical depths and timescales of radiation diffusion in the system.

To deal with this issue, we developed a four-stage method that is described below. The method first allows the disk to reach the thermal equilibrium (Sects. 2.3.1 and 2.3.2), then the planet is introduced while the disk adjusts to its presence (Sect. 2.3.3), and finally the gap is allowed to fully open and the planet-disk interactions are analysed (Sect. 2.3.4). In all stages, we simulate only one half of the disk in colatitude starting at the midplane and extending over Δϕ = 14°. We assume that the solution is symmetric with respect to the midplane, which is ensured by the reflective boundary condition for v_ϕ and the zero gradient boundary condition for the remaining quantities.

2.3.1 Hydrostatic relaxation

To find an equilibrium disk consistent with the target accretion rate Ṁ for a given constant value of α, we closely follow the hydrostatic relaxation recipe of Flock et al. (2013) while introducing slight modifications. In this section, we consider T = ɛ∕(c_Vρ) and P = (γ − 1)ɛ as independent variables (rather than ɛ and ρ).

The grid for the hydrostatic relaxation is effectively 2D (N_θ = 1) as we consider that the solution is axially symmetric in the azimuth. Moreover, the radial extent of the disk is larger than in the remaining simulation stages – we set r_min = 0.52 and r_max = 19.968 au and we resolve the radius with N_r = 374 cells. The increased radial extent improves the accuracy of the radially integrated optical depth (because in Eq. (7) we have to guess how the stellar radiation is blocked by the disk material inwards from r_min, as we explain in the following) and allows the outer disk to flare freely.

Our initial state for the hydrostatic relaxation follows the optically thin temperature (e.g. Dullemond et al. 2001), where T_thin is constant on cylinders with radius R. The density is initialized using vertically isothermal Gaussians as where z is the vertical distance from the midplane and (12)

is the target surface density given by the viscous evolution theory of razor-thin disks. The connection of the expression to our 3D model is provided through the vertically averaged viscosity (defined below). From the initial state, we iterate over the following steps.

Step I: Eqs. (5) and (6) are solved in an implicit form (Chrenko & Lambrechts 2019) while calculating the optical depth inwards from the grid as τ₀(r_min, ϕ) = κ_⋆ρ(r_min, ϕ)(r_min − 6R_⋆) (Flock et al. 2013). Since Eqs. (5) and (6) need to be advanced over a certain time step, we estimate it using a characteristic timescale of radiation diffusion Δt_dif (see Flock et al. 2013). In subsequent iterations, the time step is adaptively prolonged (or shortened) if the relative change of temperature drops below 0.1% (or exceeds 10%). At the end of Step I, new T and E_R fields are obtained and remain fixed during the subsequent steps.

Step II: We calculate a new density profile in the midplane as and convert it to P. Then we solve the equations of the hydrostatic equilibrium (e.g. Masset & Benítez-Llambay 2016)

The equations can be combined together as (15)

and solved in an implicit form by starting from the midplane pressure profile and integrating over the discrete steps Δϕ in colatitude. The implicit solution is obtained by the successive over-relaxation method with the relative precision ɛ_sor = 10⁻⁸. In this manner, a vertically stratified profile P(r, ϕ) is obtained that can be converted back to ρ(r, ϕ).

Step III: At each vertical column of cells (for a given r), we calculatethe density-weighted vertically averaged viscosity (16)

where (17)

In our calculations, the integrals are replaced with discrete sums from the midplane to the disk surface and are multiplied by a factor of two to account for both disk sides. We then recalculate Σ_target by plugging into Eq. (12) and we normalize ρ by the multiplicative factor f_norm = Σ_target∕Σ. Then the iterative procedure returns to Step I unless the relative change in T and P during a single iteration is ɛ_hst = 10⁻⁵ or smaller.

The velocity field (v_r, v_θ, v_ϕ) of a hydrostatically relaxed disk is calculated ex post. We assume that v_ϕ = 0 and obtain v_θ from Eq. (13). The remaining component v_r can be estimated from the azimuthal component of the momentum equation owing to the used assumptions of hydrostatic equilibrium (∂_t = 0), axial symmetry (∂_θ = 0), and negligible vertical motions (see Appendix A; also Takeuchi & Lin 2002; Fromang & Nelson 2006; Fromang et al. 2011; Jacquet 2013): (18)

The relevant components τ_ij of the viscous stress tensor T do not depend on v_r if the assumptionshold, thus allowing v_r to be determined. Applying Eq. (18) is important because v_r is vertically stratified in 3D disks with constant α-viscosity.

To verify that the hydrostatic stage terminated successfully, we check that the accretion rate indeed exhibits the target value Ṁ by comparing it to the instantaneous mass flux (19)

at all radii.

2.3.2 Hydrodynamic relaxation

We perform a hydrodynamic relaxation of the disk as the second stage of our simulations. The goal is to verify that (i) the disk structure does not significantly change; (ii) the accretion rate remains close to Ṁ when full viscous stresses are introduced. Starting from the final state of the hydrostatic relaxation, we truncate the disk to the radial extent given in Table 1 but we still maintain the assumption of axial symmetry. We then solve the full set of Eqs. (3)–(6) over the timescale of 6000 P_orb where P_orb is the orbital period at 5.2 au.

Special boundary conditions are adopted for independent variables ρ, ɛ, E_R, and v. At radial boundaries, we use a Keplerian extrapolation for v_θ (e.g. Bitsch et al. 2014) and the zero gradient condition for E_R, v_ϕ, and v_r. The latter is reflected if an inflow into the domain is detected with a Mach number 0.1 or larger (Flock et al. 2013). Using a radial power-law extrapolation for T, we calculate (Eq. (16)) as well as Σ_target (Eq. (12)) in each ghost ring. We copy ρ from the first active cell, calculate Σ (Eq. (17)), and rescale ρ by f_norm = Σ_target∕Σ. The rescaling ensures that the boundary accretion rate of the disk remains consistent throughout the simulation. Finally, we compute ɛ = ρc_VT in each ghost cell.

Additionally, the boundary conditions are supplemented with the wave-killing zones of de Val-Borro et al. (2006) in the radial direction and near the disk surface (not near the midplane). We damp only the velocity components in a way that v_θ is damped towards its azimuthal average while v_r and v_ϕ are damped towards their hydrostatic values. After several initial tests, we chose a rather stringent damping timescale equal to 0.03 of the local Keplerian period.

The described boundary conditions are also used in the remaining simulation stages and so is τ₀ (r_min), which we prescribe according to the value extracted from the hydrostatic stage (at the 20th radial ring of the hydrostatic grid).

2.3.3 Planet insertion

Next, the grid is expanded in the azimuth to the final number of N_θ = 628 zones. The arrays of quantities obtained during the hydrodynamic relaxation are copied azimuthally to cover the expanded grid. The disk is then evolved for t = 0–500 P_orb (here we define the time origin t = 0) during which we insert the planet into the simulation. The planet mass is gradually increased from zero to its final value during t = 0–250 P_orb.

As the planet grows, it starts to open the gap. The gas is pushed away from the planetary orbit and then it continues to spread viscously from the gap edges. The process tends to be slow and it is beneficial to speed it up in numerical simulations by allowing the planet to accrete gas (Crida & Bitsch 2017). The accretion is achieved by removing the fraction of gas Kf(d)δt from within the Hill sphere (Kley 1999) where δt is the hydrodynamic time step determined by the Courant–Friedrics–Lewy (CFL) condition of FARGO3D, K is an arbitrary parametrization of the accretion efficiency, and (Crida et al. 2016) (20)

We let K to increase from 0 to 1 over 250 P_orb and then we decrease it back to 0 during t = 250–500 P_orb. When K > 0, the planet is assumed to behave as a mass sink; the gas is removed from the simulation but it is not added to the planet mass, nor is the gas momentum. Apart from the planet insertion, the planet is typically non-accreting (K = 0), unless stated otherwise.

2.3.4 Main stage

After the planet insertion, the simulation is continued until the measured disk-driven torque converges to a stationary value. The convergence is only achieved once the gap profile is settled and the disk structure becomes adjusted to its presence. The planet is kept on a fixed circular orbit and the obtained torque is therefore the static torque.

The main stage typically covers t = 500–3500 P_orb. The calculation is numerically demanding because (i) the 3D grid has a relatively large number of cells (≃ 5.4 × 10⁶); (ii) fast wave propagation in low-density regions of the disk diminishes the maximum allowed time step through the CFL condition. The latter becomes especially restrictive once the gap is opened and large density contrasts are produced between the midplane and the disksurface. This is partially compensated for by the wave-killing procedure and by introducing a volume density floor ρ_floor = 10⁻²² g cm⁻³.

Our simulations were run on CPU clusters NASA Pleiades and IT4I Salomon. To achieve a reasonable speedup, weused a domain decomposition and a hybrid parallelization based on the Message Passing Interface (MPI) and Open Multi-Processing (OpenMP). A single simulation was usually spawned over ≃ 550 CPU cores and consumed ~ 10⁵ CPU hours.

2.4 Reference non-radiative simulations

Since we aim to isolate the influence of gap irradiation on the Type II torque, it is beneficial to compare the results of simulations with stellar irradiation (Sect. 2.3) to reference simulations that preserve the equilibrium disk temperature even after the gap opening (i.e. they neglect the increased amount of stellar heating at the outer gap edge). For our reference model, we replace Eqs. (5) and (6) with a single energy equation that neglects any radiative effects: (21)

where the subscript “0” stands for quantities at t = 0 (at the beginning of planet insertion) and t_cool is the cooling timescale, which is set to 10⁻³ of the local orbital period. Reference simulations begin with the planet insertion stage (as there is no need to recalculate the unperturbed disk). Due to the short cooling timescale, the reference model is expected to behave similarly to 3D locally isothermal simulations.

Fig. 1 Radial profiles of the accretion mass flux FM (top), aspectratio h (middle), and surface density Σ (bottom) after the hydrostatic relaxation (dashed blue curve) and hydrodynamic relaxation (solid black curve). The latter spanned 6000 Porb. The red dotted curve shows a fit of the h ∝ r2∕7 dependence to the result of the hydrodynamic relaxation. The profiles demonstrate that our equilibrium disk has its global thermodynamics governed by passive heating (stellar irradiation) and its accretion rate is almost uniform.

3 Results

3.1 Fiducial case

Here, we analyse the simulation based on our fiducial parameters. We study the disk structure (Sect. 3.1.1), planet-induced perturbations (Sects. 3.1.2 and 3.1.3), torque measurements (Sect. 3.1.4), and gap edge instabilities (Sects. 3.1.5 and 3.1.6).

3.1.1 Equilibrium disk

Figure 1 shows the equilibrium radial profiles of F_M, h, and Σ at the end of the hydrostatic and hydrodynamic relaxation (Sects. 2.3.1 and 2.3.2, respectively). After the hydrostatic relaxation, the radial mass flux due to disk accretion is very close to the target value of Ṁ = 10⁻⁸ M_⊙ yr⁻¹. During the hydrodynamic relaxation, the introduction of the full viscous stress results in small structural changes with respect to the hydrostatic state. Since the radial gas velocity v_r is sensitive even to small perturbations (e.g. Dürmann & Kley 2015), the mass flux departs from the target valueand then slowly converges to the state depicted in Fig. 1. Although the final F_M (r) does not perfectly match the hydrostatic state, we consider the differences acceptable because no significant departures are apparent in h(r) and Σ(r).

The displayed profile of h(r) is computedfrom the midplane temperature as (e.g. Bitsch et al. 2014) (22)

Clearly, h(r) corresponds to a passively heated protoplanetary disk because it can be well characterized by a least-squares fit of h ∝ r^2∕7 in accordance with Chiang & Goldreich (1997). The Σ(r) profile exhibits Σ_p = 385 g cm⁻² at the planet location and the characteristic disk mass (Dürmann & Kley 2015) (23)

attains M_D = 1.17 M_J for the fiducial set of parameters.

Fig. 2 Temperature T0 before planet insertion (top) and the temperature difference T − T0 after gap opening (bottom; taken at t = 3500 Porb) in the vertical plane of the disk. The quantities are azimuthally averaged. The Hill sphere and the location of a Jupiter-mass planet are displayed in the bottom panel. We point out that the scale of the vertical axis is exaggerated. The temperature increase due to the irradiation of the outer gap edge is clearly visible in the bottom panel.

3.1.2 Temperature perturbation due to gap edge irradiation

Once the planet is inserted into the disk, it starts to open the gap. Due to gas clearing in the gap centre, the outer gap edge becomes more exposed to stellar irradiation and a large-scale temperature variation is expected. Figure 2 shows the azimuthally averaged temperature distribution in the meridional plane. Before planet insertion (top panel of Fig. 2), the temperature map exhibits features typical for passive protoplanetary disks (e.g. Flock et al. 2013, 2017). Two distinct disk layers can be distinguished – a hotter photosphere and a cooler interior. The temperature rise in the photosphere appears because the region is optically thin to stellar irradiation (τ_⋆ ≲ 1). The interior below the irradiated surface, on the other hand, is only heated by the reprocessing of the thermal radiation and becomes nearly vertically isothermal.

By studying the temperature variation after the gap opening (bottom panel of Fig. 2), one can see thatthe embedded planet substantially changes the thermal structure of the surrounding disk. An overheated layer appears above the protoplanet since both the gap clearing and vertical disk contraction increase the extent of the photosphere. Yet another temperature excess appears as a hot column spanning r ≃ 6–8 au ≃ 1.2–2 r_p. This excess is due to the increased amount of irradiation intercepted by the exposed outer gap edge. The minor temperature excess at around ≃ 2 au is a leftover from the transitional phase during which the gas is expelled away from the planetary orbit (the excess slowly disappears over time).

Conversely, there is a temperature deficit in the lower layers of the disk outwards from the inner gap edge and also in the upper layersoutwards from the puffed up outer gap edge. In both cases, the respective region is shielded from direct stellar illumination and remains heated only by radiative diffusion. Overall, the thermal structure after gap openingis in a good agreement with the results of Jang-Condell & Turner (2013) who reported the same features.

3.1.3 Disk structure after gap opening

The response of the global disk structure to the gap opening and temperature variations is shown in Fig. 3. The perturbed surface density Σ∕Σ_p of the disk reveals that, for the given combination of parameters, the gap becomes relatively deep (compare e.g. to Dürmann & Kley 2015; Robert et al. 2018). To assess if there are any peculiarities in the gap profile, we compared the gap width and depth with predictions resulting from the 2D locally isothermal simulations of Kanagawa et al. (2016). They derived a gap width of (24)

and depth of (25)

where (26)

By plugging in the values from our simulation, we find the gap depth to be Σ_min∕Σ₀ ≃ 1.2 × 10⁻³, which corresponds very well to the density drop in the gap centre (as indicated by the dotted horizontal line in Fig. 3). Regarding the gap width, Kanagawa et al. (2016) define it as the radial extent where the azimuthally averaged surface density is smaller than half of the initial surface density. In our simulation (see the hatched band in Fig. 3), the definition² leads to Δ_gap ≃ 0.78 r_p, which is best recovered if the constant of proportionality in Eq. (24) is 0.35, only slightly smaller than 0.41 derived by Kanagawa et al. (2016). Therefore, the gap opened in our 3D radiative disk is similar to a 2D situation, as already pointed out by Fung & Chiang (2016). The only notable feature of the gap is a slight asymmetry – the inner half of the gap is more depleted compared to the outer half.

Turning our attention to the aspect ratio h (middle panel of Fig. 3), the most prominent feature of the final state with respect to t = 0 is the bump that peaks close to the outer gap edge. Since h (and c_s) scales with (Eqs. (2) and (22)), we can deduce that any variations of h reflect the perturbed thermal structure (Fig. 2). Specifically, the bump at the outer gap edge arises because the local heating becomes more efficient after the gap opening. The shade of the inner gap edge is responsible for the drop of h at r ≃ 0.75 r_p and the puffed-up outer gap edge shadows the region at r ≳ 2.5 r_p.

The profile of the midplane sound speed perturbation (bottom panel of Fig. 3) provides a useful comparison to Hallam & Paardekooper (2018). Since they used a 2D vertically averaged model without radiation physics, they had to estimate the increase in the sound speed due to edge illumination. By treating the outer gap edge as a disk rim, they deduced a boost of c_s by a factor of 1.8–2.5 in their fiducial case. The authors pointed out that this is rather an upper limit and that a realistic boost would likely be less strong due to the blocking of the starlight by the inner disk. Indeed, our fiducial simulation reveals that the boost of c_s assumed by Hallam & Paardekooper (2018) was probably an overestimate since we measure the maximum increase as being by a factor of ≃ 1.2. It would be difficult to achieve a factor of 2 because T would have to rise by a factor of 4, from ≃ 50 K to at least ≃ 200 K, which is hotter than the photosphere at the location of the outer gap edge in our fiducial case.

Finally, since Hallam & Paardekooper (2018) used a Gaussian function to mimic the sound speed perturbation, it is worthwhile checking how well the Gaussian represents the peak of the c_s ∕c_s,0 profile. As shown in Fig. 3, a good representation for r > r_p can be obtained by a least-squares fit of the skewed Gaussian (27)

leading to the tail value , amplitude , central position r_G = 1.23 r_p, standard deviation σ = 0.56, and skewness ζ = 4.65. By setting and ζ = 0, one recovers the exact form of the Gaussian assumed by Hallam & Paardekooper (2018) and the least-squares fitting of the remaining free parameters then leads to , r_G = 1.55 r_p, and σ = 0.24, implying the full width at half maximum . But such a Gaussian cannot properly reproduce the skewed shape nor the tails of the peak in Fig. 3.

Fig. 3 Azimuthally averaged profiles of the surface density perturbation Σ∕Σp (top), aspectratio h (middle), and midplane sound speed perturbation cs∕cs,0 (bottom) during the course of the fiducial simulation with irradiation. Solid curves corresponding to various simulation times are distinguished by colour. The vertical dashed line is the planet location, the grey vertical band marks the Hill sphere, and the hatched vertical band shows the extent of the gap width (following the definition of Kanagawa et al. 2016). The horizontal dotted line marks the gap depth according to Kanagawa et al. (2016). Dashed curves in the bottom panel are the least-squares fits of the Gaussian function of Hallam & Paardekooper (2018; magenta) and the skewed Gaussian (turquoise) given by Eq. (27). The peak of h and cs can be directly related to the increased heating of the outer gap edge by stellar irradiation (see Fig. 2).

Fig. 4 Temporal evolution of the normalized gas-driven static torque γΓ∕Γ0 in the fiducial case. We display the unfiltered measurement with the output sampling of 1∕20 Porb (grey curves), as well as the moving average of the one-sided inner, outer, and total torque (blue, red, and black curve, respectively). The horizontal green dashed line shows the value of the converged torque (averaged over the last 1000 Porb). The final value is γΓ∕Γ0 = −0.035.

3.1.4 Torque evolution

Figure 4 shows the temporal evolution of the static torque exerted by the disk on the planet. We normalize the torque Γ as γΓ∕Γ₀ where (e.g. Paardekooper et al. 2010) (28)

and all quantities correspond to the state before planet insertion. The quantity Γ₀ is defined to reflect the basic dependencies (e.g. ; ~ h⁻²) of Type I torques but, for the sake of consistency, it is usually used to characterize Type II migration as well (see Kanagawa et al. 2018). When normalizing the torque as γΓ∕Γ₀, the factor γ accounts forthe difference in the sound speed between isothermal and non-isothermal models (Baruteau & Masset 2008).

From Fig. 4, one can see that the total torque undergoes fast oscillations related to the outer one-sided torque. These oscillations first appear during planet introduction (t <500 P_orb) simply because the disk undergoes abrupt changes as the gap is being opened. During t ≃ 700–2000 P_orb, oscillations appear again in episodes with a relatively small and gradually decreasing amplitude. At about t ≃ 2000 P_orb, the amplitude of torque oscillations increases and they no longer vanish. The origin of torque oscillations during the main simulation stage (t > 500 P_orb) is not clear at first glance and will be investigated later.

To filter out fast torque oscillations, we smoothed out the time series of our torque measurement by a moving average with a window size of 50 P_orb. To gain the final value of the torque, we calculated the arithmetic mean over the last 1000 P_orb of our simulation and we also verified that prolonging the simulation to 5000 P_orb does not lead to a substantially different total torque. We measured γΓ∕Γ₀ = −0.035.

The question now arises – is the measured torque reduced due to gap edge irradiation? To answer the question, we performed a reference non-radiative simulation that neglects gap illumination (Sect. 2.4) and preserves the unperturbed thermal structure of the disk (as in the top panel of Fig. 2). We obtained the converged total torque which implies that for our fiducial parameters, gap irradiation does not reduce the magnitude of the total torque³.

To provide further insight into where the torque is generated, in Fig. 5 we plot the radial distribution of the specific (per unit mass) torque Γ(r) and the cumulative torque Γ(< r). The former basically measures the torque exerted on the planet by the gas located on a grid annulus with the radius r; the latter represents the total torque summed from r_min to r. The profile of Γ(r) for the simulation with irradiation is obtained by calculating the time average over 10 P_orb.

The profile of Γ(r) (top panel of Fig. 5) reveals that the largest difference in the presence of gap irradiation appears for several peaks just outside the outer gap edge. However, this difference has a negligible influence on the total torque because the peaks have a tendency to average out, as apparent from the Γ(< r) profile (bottom panel of Fig. 5) – the change in the cumulative torque between r ≃ 1.4 r_p and r ≃ 3 r_p is rather small. The greatest gain in the cumulative torque appears across the inner half of the gap, and the greatest loss appears across the outer half. But in this region, the profiles of Γ(r) measured with and without gap irradiation are qualitatively very similar.

So far, we can see that the obtained result is altogether negative; there seems to be no significant influence of gap irradiation on the total torque. But we will demonstrate later in Sect. 3.2 that the total torque can be reduced when planetary ordisk parameters are changed.

Finally, let us point out that Γ(r) in Fig. 5 shows that torque oscillations arising in our fiducial simulation originate in the outer half of the gap. The rest of Sect. 3.1 is dedicated to finding the source of the torque oscillations.

Fig. 5 Radial profile of the specific torque Γ(r) (top) and the cumulative torque Γ(<r) (bottom) for the fiducial set of parameters. The simulation with irradiation (solid black curve) and the reference non-radiative simulation (dotted blue curve) are shown. The grey curves represent 200 samples of Γ(r) recorded over 10 Porb that were used to calculate the black curve as the arithmetic mean and they trace the torque oscillations. The planet location andgap width are shown as in Fig. 3.

Fig. 6 Temperature difference in the disk midplane with respect to the unperturbed state (similar to the bottom panel of Fig. 2). The figure is taken from our fiducial simulation at t = 3500 Porb and revealsstreamers (filaments) stretching across the gap and spiral arms, which are offset with respect to the planet-induced perturbations. Some of these structures are marked with arrows.

3.1.5 Outer gap edge instability

Torque oscillations can only be induced by azimuthal asymmetries of the gas distribution that are not corotating with the planet. In the case of giant planets, such asymmetries are usually caused by the excitation of vortices either in the co-orbital region or at gap edges (Koller et al. 2003; Li et al. 2005; de Val-Borro et al. 2007; Ou et al. 2007; Lin & Papaloizou 2010; Les & Lin 2015). Indeed, Fig. 6 shows the temperature perturbation in the disk midplane and contains overheated streamers (filaments) (Fung & Chiang 2016) and additional spiral wakes excited mostly in the outer half of the gap. It is natural to assume that these structures trace the presence of vortices.

Figure 7 compares the perturbed surface density distribution between the reference simulation and the irradiated simulation with fiducial parameters. In the latter case, the outer gap edge is clearly perturbed in a wave-like manner and one of the additional spiral wakes centred at r ≃ 1.2 r_p, θ ≃ −1 rad is visible. We identified that the wavy perturbation of the outer gap edge is not static in the reference frame of the planet and thus it can only be caused by vortical structures propagating at a non-zero phase speed.

The final demonstration of the presence of vortices is provided in Fig. 8 where we study the relative perturbation of the vorticity component (29)

calculated from the 2D midplane velocity field in the frame corotating with the planet. In Fig. 8, we see that the vorticity is distributed in an orderly fashion in the reference simulation where there are neighbouring sheets of large positive and negative vorticity that delimit the co-orbital region. Such a vorticity distribution no longer exists in the simulation with stellar irradiation where it is disturbed by the presence of several vortices.

Before proceeding with the (in)stability analysis, let us point out that Figs. 6 and 7 contain additional valuable information. Specifically, Fig. 6 reveals azimuthal asymmetries in the temperature variations from which we calculated the maximum temperature contrast reached in the planetary spiral wake with respect to the disk background as 17%. The obtained value is in a good agreement with Ziampras et al. (2020; their contrast from 2D simulations was 15%). Figure 7 then shows that in our radiative simulation, the spiral arms have a decreased density contrast (again as in Ziampras et al. 2020) and their winding is less tight in the outer disk. The latter can be explained by the dependence of the pitch angle (Zhu et al. 2015), which grows as h becomes puffed up in the outer disk in the presence of gap irradiation. The differences in the planetary wake are responsible for the differences found in Fig. 5.

Fig. 7 Surface density perturbation (Σ − Σ0(r))∕Σ0(r) relative to the equilibrium state before planet insertion. The figure is taken at t = 3500 Porb and shows the reference non-radiative simulation (top) and the simulation with stellar irradiation (bottom) for the fiducial set of parameters. The outer gap edge in the bottom panel exhibits additional azimuthal asymmetries.

Fig. 8 Midplane vorticity perturbation (ω⊥− ω⊥,0)∕ω⊥,0 relative to the equilibrium state. The figure is taken at t = 3500 Porb and shows the reference non-radiative simulation (top) and the simulation with stellar irradiation (bottom) for the fiducial set of parameters.

3.1.6 Stability analysis

The excitation of vortices that ultimately lead to torque oscillations has to arise due to a hydrodynamic instability. In this section, we analyse the vulnerability of the disk to the most common hydrodynamic instabilities that can occur in the presence of embedded planets: the Rayleigh instability, the buoyant instability, and the Rossby wave instability.

The Rayleigh instability can occur at gap edges where planet-induced perturbations of the pressure significantly modify the local orbital velocity (Kanagawa et al. 2015; Fung & Chiang 2016). In stable disks, the angular momentum per unit mass j = R²Ω increases with radius (Chandrasekhar 1961). For the Rayleigh stability criterion, we use (30)

where the brackets denote the vertical density-weighted average (as in Eq. (16)).

The buoyant instability appears in non-barotropic disks wherever there is a misalignment between pressure and density gradients (e.g. Klahr & Bodenheimer 2003; Petersen et al. 2007; Lesur & Papaloizou 2010), which is often satisfied in planet-driven shocks (e.g. Ou et al. 2007; Richert et al. 2015). To assess the buoyant stability, we adopt a form of the Solberg-Høiland criterion (Rüdiger et al. 2002) (31)

where we calculate the square of the epicyclic frequency as (32)

and the square of the radial Brunt-Väisälä frequency as (33)

where Π is the vertically integrated pressure and γ_2D = (3γ − 1)∕(γ + 1) is the 2D adiabatic index (Klahr 2004).

The Rossby wave instability (Lovelace et al. 1999; Li et al. 2000) is a result of the velocity shear at the edges of steep planet-induced gaps. It can become excited when the generalized potential vorticity (Lin 2013) (34)

develops a local inflection point. The quantity ξ essentially represents an entropy-modified version of vortensity (vorticity divided by surface density). Inflection points of ξ can arise (Koller et al. 2003; Lin & Papaloizou 2010) because the disk material crossing the shock associated with the planetary spiral wake undergoes a modification of the vortensity (Li et al. 2005) as well as of the entropy (because there is a temperature jump across the shock; see Fig. 6 and Les & Lin 2015). Typically, two coupled Rossby waves are excited around an inflection point of ξ and they emit spiral density waves (Meheut et al. 2010).

The stability analysis is summarized in Fig. 9 and implies that the disk remains stable to the Rayleigh and buoyant instabilities because the respective criteria are positive at all radii. However, our fiducial simulation with gap irradiation exhibits an inflection point of ξ at r≃ 1.17 r_p and is therefore susceptible to the Rossby wave instability and we identify it as the source of vortices. Using the Hill sphere width to guide the eye, we can conclude that the inflection point appears at the boundary between the librating and horseshoe streamlines (because the width of the horseshoe region is ≃ 2.45 R_H for giant planets; Masset et al. 2006a). In the absence of gap irradiation (as investigated by our reference simulation), no vortices are excited because there is no inflectionpoint of ξ.

Although the Rossby wave instability is a viable explanation, we stress that our analysis might not be entirely conclusive. For example, the Rossby wave instability has often been found to produce a single merged vortex (e.g. Les & Lin 2015), which we do not see in our simulations. Additionally, we detect vortices even for viscosity values for which they were previously found to dissipate (Fu et al. 2014). It is thus possible that our grid resolution is not sufficient to properly capture the behaviour of vortices. But it is also possible that we see realistic effects related to the behaviour of vortices in 3D or to the influence of gap irradiation. Further investigation is beyond the scope of this paper.

Fig. 9 Azimuthally averaged radial profiles of various hydrodynamic stability criteria in the vicinity of the planet. From top to bottom, we show the Rayleigh criterion (the left-hand side of Eq. (30)), Solberg-Høiland criterion (the left-hand side of Eq. (31)), and the generalized potential vorticity (Eq. (34)). The simulation with gap irradiation (solid black curve) and the reference non-radiative simulation (dotted blue curve) are compared. The grey rectangle shows the extent of the Hill sphere. The inflection point of the black curve at about r ≃ 1.17 rp in the bottom panel suggests that the disk is Rossby-unstable.

3.2 Dependence on parameters

In this section, we perform a coarse parametric study to test the behaviour of the Type II torque in stellar-irradiated disks under various conditions. We explore the dependence on the planet mass M_p, α viscosity, and disk accretion rate Ṁ. Since it would be numerically expensive to sample mutual combinations of these parameters in 3D, we always vary a single parameter at a time while keeping the others fixed to their fiducial values. The only exception are simulations with varying α. For a fixed value of Ṁ, a variation of α would change Σ of the disk (Eq. (12)) and thus also M_D (Eq. (23)). But both the disk mass and viscosity can affect Type II migration (Dürmann & Kley 2015; Robert et al. 2018) and thus it is undesirable to mix these two effects together. To avoid this, we adjust Ṁ when varying α in order to keep M_D fixed. We point out that the dependence on the disk mass is studied separately in the simulation set with variable Ṁ (and other parameters fixed).

The summary of simulations performed in our parametric study is given in Table 2. The simulation time span was usually t = 3500 P_orb, with the exception of simulations M_p = 0.5 M_J and Ṁ = 10⁻⁹ M_⊙ yr⁻¹, which were prolonged to t = 5000 P_orb to improve the convergence of the torque. For cases M_p = 0.1, 0.18, 0.25, 0.5, 2 M_J and α = 5 × 10⁻⁴, 5 × 10⁻³, we also performed reference simulations without gap irradiation.

3.2.1 Summary of the static torque measurements

Figure 10 summarizes the main results of our parametric study (including the fiducial simulation), namely the dependence of the static torque on M_p and α. Reference simulations that neglect gap edge irradiation are shown as well. In accordance with previous studies (e.g. Dürmann & Kley 2015), the static Type II torque measured in reference simulations is always negative and its normalized magnitude grows with (i) decreasing planet mass and (ii) increasing α viscosity.

When gap irradiation is taken into account, Fig. 10 reveals that the magnitude of the torque is reduced for M_p < 1 M_J and α >10⁻³. In the case of M_p = 0.25 M_J, the torque even switches its sign from negative to positive, which would result in outward migration. To quantify the torque reduction, we compute the reduction factor f_irr in (35)

where Γ_ref is the torque measured in reference non-radiative simulations. Values of f_irr (if available) are displayed in Fig. 10 as labels next to respective data points. Overall, Fig. 10 constrains the parametric space in which the mechanism of Hallam & Paardekooper (2018) becomes important.

The behaviour of trends in Fig. 10 suggests that the efficiency of the torque reduction can be related to the gap depth and width. The latter is captured in Fig. 10 by the additional horizontal axis, which shows the position of the outer gap edge r_edge. As a reminder, let us recall that planetary gaps become narrower and shallower with (i) decreasing planet mass or (ii) increasing viscosity (in both cases, the planet gets less efficient in expelling the gas away from its orbit).

For cases M_p ≥ 1 M_p and α ≤ 10⁻³, which correspond to the widest and deepest gaps obtained in our study, the trends behave similarly to the reference simulations. However, the trends exhibit turnover points at M_p = 1 M_J and α = 10⁻³ (which exactly corresponds to our fiducial case) and their dependence on the respective parameter becomes reversed with respect to the reference curves for 0.25 M_J < M_p < 1 M_J and α > 10⁻³. In this range of parameters, the torque reduction clearly becomes more efficient for shallower and narrower gaps. For M_p < 0.25, the dependence on the planet mass switches back to the expected trend and the effect of the torque reduction starts to weaken (the “black” trend starts to bend back to the “grey” trend).

Based on the γΓ∕Γ₀(M_p) dependence, we can summarize that the torque felt by a giant planet is reduced due to stellar irradiation for moderately wide and deep gaps. The torque reduction vanishes once the gap becomes too wide and deep (as studied in Sect. 3.2.2) or too narrow and shallow (because the temperature excess outwards from the planet disappears if M_p becomes too small). The γΓ∕Γ₀(α) dependence exhibits the same behaviour but we can only see the onset of the torque reduction for moderately wide gaps. In order to recover the cutoff observed for narrow gaps, we would have to test even larger values of α but it is reasonable to assume that the “black” trend would once again bend back towards the “grey” trend.

Additionally, Fig. 10 seems to support our claim that torque oscillations do not affect the mean value of the torque since there are no apparent changes in the displayed trends at the transitions between simulations with and without instabilities (as distinguished by filled and open symbols, respectively).

The dependence of the Type II torque on Ṁ is not shown in Fig. 10 because of numerical issues that we encountered for the extremal values of the accretion rate and which would most likely lead to an erroneous torque value. For Ṁ = 10⁻⁹ M_⊙ yr⁻¹, the planet is considerably more massive than the disk (M_D = 0.12 M_J) and consequently, the planet-induced perturbations of the outer disk are so strong that neither the disk structure nor the torque can converge before 5000 P_orb. For Ṁ = 10⁻⁷ M_⊙ yr⁻¹, the disk is relatively massive (M_D = 10.7 M_J) and its optically thick interior is vertically too expanded and the disk photosphere too thin with respect to the opening angle of our computational domain.

However, the remaining two simulations with Ṁ = 5 × 10⁻⁹ and 5 × 10⁻⁸ M_⊙ yr⁻¹ proceeded normally and we obtained the normalized value of the torque as γΓ∕Γ₀ = −0.032 and − 0.037, respectively.Both these values are very similar to the fiducial simulation with γΓ∕Γ₀ = −0.035, which suggests that the measured torque has a very weak dependence on the disk mass. But we admit that such a statement is based on a very small sample of simulations and requires future investigation.

Fig. 10 Dependence of the normalized torque γΓ∕Γ0 on the planet mass Mp (top) and α viscosity (bottom). Additionally, the top horizontal axis of each panel shows the position of the outer gap edge redge as resulting from Eq. (24). Circles and solid black lines correspond to simulations with gap irradiation; diamonds and dashed grey lines correspond to reference non-radiative simulations. Open symbols correspond to simulations without torque oscillations (i.e. Rossby-stable) while filled symbols mark simulations with gap edge instabilities. Data points are labelled with the corresponding torque reduction factor firr (Eq. (35)). Clearly, gap edge irradiation reduces the magnitude of the torque (or even causes torque reversal) for Mp < 1 MJ and α > 10−3.

Fig. 11 Radial profiles of the surface density perturbation Σ∕Σp (top), temperature T (middle), and specific torque Γ(r) (bottom). Simulations with different planet masses Mp = 0.25 (red), 0.5 (blue), and 1 MJ (black) are shown. Simulations with irradiation and reference non-radiative simulations are distinguished using solid and dotted lines,respectively. The top two panels are azimuthally averaged.

3.2.2 Torque reduction with decreasing planet mass

Let us investigate how the torque reduction operates in the mass interval between M_p = 0.25 (for which the torque–mass dependence in Fig. 10 peaks) and the fiducial case. In this interval, the torque-mass dependence exhibits an unexpected ‘reversed’ trend and it is thus worth a deeper analysis. We focus on cases M_p = 0.25 and 0.5 M_J and we compare them to the fiducial case M_p = 1 M_J. Figure 11 shows the resulting radial profiles of the surface density perturbation Σ∕Σ_p, temperature T, and specific torque Γ. In all cases, the modification of the disk structure is qualitatively similar. As the planet mass decreases, the gap becomes narrower and shallower as one would expect (Eqs. (24)–(26)). Consequently, the peak of the temperature excess behind the outer gap edge decreases and also recedes inwards, closer to the planet.

The torque distribution exhibits the most prominent changes in the gap region as well. With decreasing planet mass, the amplitude of the peaksclosest to the planet increases as a greater amount of gas remains in the gap region and the peaks themselves shift towards the planet. As there are no qualitative differences that one could easily recognize in the shape of individual curves of Γ(r) in Fig. 11, the torque reduction identified in Fig. 10 can only be explained by the asymmetry of the one-sided torques.

As a confirmation, we compare the one-sided torques in Fig. 12 where we show Γ(r) of the inner and outer peak closest to the planet as a function of the distance from the planet. Additionally, we normalize the torque so that the inner peak always has a maximum of Γ = 1. With this normalization, the torque reduction becomes apparent. In the reference non-radiative simulations, the inner-outer asymmetry is always such that the outer negative one-sided torque dominates. However, in simulations with gap irradiation and M_p = 0.25 and 0.5 M_J, the extremal values of the one-sided torques are almost identical. For M_p = 0.5 M_J, the outer one-sided torque has a somewhat wider profile than the inner one and the total torque thus remains negative, albeit with a reduced magnitude. For M_p = 0.25 M_J, the one-sided torques in the depicted interval of Δr almost cancel out. Interestingly, we notice that the inner-outer asymmetry is reduced even for the fiducial case (which was not apparent from Fig. 5) although the effect is not sufficient to substantially change the total torque.

Finally, Fig. 13 shows the cumulative specific torque for the discussed cases. Clearly, the largest variation of Γ(r) arises approximately from the region 0.65 r_p < r < 1.35 r_p, which makes the range of Fig. 12 justifiable. Figure 13 reveals that without gap irradiation, the cumulative torque across the outer half of the gap drops to almost the same specific value regardless of the planet mass. Once gap irradiation is considered, the drop across the outer half of the gap becomes less strong with decreasing planet mass.

Fig. 12 Specific torque Γ(r) as a function of the radial separation from the planet for various planetary masses (as in Fig. 11). The predominantly positive branch of the curves corresponds to the inner one-sided torque (arising from the interval 0.65 rp < r < 1 rp) and the negative branch is the outer one-sided torque (arising from the interval 1 rp < r< 1.35 rp). The aim of the figure is to assess the inner-outer asymmetry of the one-sided torques and Γ(r) is therefore artificially normalized to the maximum value of the inner peak. The figure reveals how gap irradiation reduces the outer Lindblad torque with respect to the inner one.

Fig. 13 Specific cumulative torque Γ(< r) as a function of the radial distance for various planetary masses (as in Fig. 11). Simulations with gap irradiation (top) and reference non-radiative simulations (bottom) are compared. The figure demonstrates how gap irradiation modifies the drop of the cumulative torque across the outer half of the gap.

Fig. 14 As in Fig. 12 but for simulations with different values of the α viscosity 10−3 (black), 3 × 10−3 (blue), and 5 × 10−3 (red).

3.2.3 Torque reduction with increasing α viscosity

For the sake of completeness, here we provide details of the torque reduction with increasing α viscosity. Figure 14 is analogous to Fig. 12 and compares the inner-outer asymmetry of the one-sided torques for the fiducial case (α = 10⁻³) and cases α = 3 × 10⁻³ and 5 × 10⁻³. Let us focus on the latter case for which the torque reduction is the strongest (see Fig. 10) and for which we also have a reference simulation to compare it to. For α = 5 × 10⁻³, we identify that (i) the inner-outer asymmetry of the extremal values is the same with and without gap irradiation; (ii) the peak of the outer torque is narrower when gap irradiation is considered (there is a wider separation between the solid and dashed red line for the tail of the outer torque than for the tail of the inner torque). Although the features identified here are more subtle compared to Sect. 3.2.2, they again demonstrate the reduction of the inner-outer asymmetry of the one-sided torques. Here the asymmetry is smeared out with increasing α viscosity.

4 Discussion

4.1 Caveats and future work

According to our analysis, the torque modification due to gap edge irradiation critically depends on the amplitude and location of the temperature perturbation related to the gap shape (as already suggested in Hallam & Paardekooper 2018). It is therefore natural to ask whether our conclusions can be generalized when (i) the planet is placed elsewhere, or (ii) a different disk model is considered, or (iii) the planet is allowed to migrate freely. Here we discuss the limitations of our model and speculate about possible implications. Their verification is left for future work.

Regarding (i), for a fixed planet mass and our fiducial disk model, the torque reduction might actually vary with the semi-major axis of the planet. Since the aspect ratio of a passively heated disk grows as h ∝ r^2∕7, the relative gap width Δ_gap∕r_p (Eq. (24)) and depth Σ_min∕Σ₀ (Eq. (25)) decreases with increasing r. Therefore, our results cannot be directly scaled to an arbitrary planetary semi-major axis and some trends can be expected. For example, the torque acting on a Jupiter-mass planet in our fiducial setup is not modified due to gap edge irradiation. But placing the planet at larger r_p where h is also larger, the gap would become narrower and shallower (measured relatively to r_p and Σ₀) and the temperature excess would be shifted as well. The Jupiter-mass planet could become affected by the torque reduction in such a situation.

Concerning (ii), our model of a passively heated irradiated disk neglects any viscous heating. As shown in Appendix B, viscous heating would make the disk hotter and h(r) would generally increase, mainly in the inner disk (within several au) where viscous heating dominates over stellar irradiation (e.g. Bitsch et al. 2014). The region around r_p = 5.2 au would exhibit a moderate boost of h(r). Consequently, a given planetary mass would open a shallower and narrower gap, with similar implications as discussed above.

The thermal balance of the disk also depends on the opacity. Combining viscous heating with our fiducial opacity model, h(r) would remain a monotonically increasing function of radial distance, albeit less steep than with stellar irradiation only. The disk surface would be irradiated under a different grazing angle and thus the temperature excess after gap opening could be slightly shifted. For an opacity law with opacity transitions (such as the frequently used law of Bell & Lin 1994), a bump of h(r) would most likely exist near the water evaporation line (Bitsch et al. 2014) and it would shield the adjacent outer part of the disk (up to ≃ 8 au) from direct stellar illumination, which might possibly prevent gap irradiation for planets located there. It remains to be studied how the torque reduction operates near such self-shadowed disk regions. Similarly, it is necessary to explore if viscous heating itself can create considerable temperature perturbations at the gap edges that would affect the inner-outer asymmetry of the one-sided Lindblad torques.

As for (iii), the static torque measured in our simulations should not be directly interpreted as an accurate representation of the migration rate. This has been pointed out by many recent studies (see Sect. 1). Once the planet is released, it becomes offset with respect to the gap centre, depending on the torque that it felt initially, and a new torque balance is established. At the same time, the accreting disk adapts to the movement of the planet, gap-crossing flows (if present) reorganize compared tothe static situation, and the gap reshapes with a certain lag. The migration then becomes coupled to the accretion flow of the disk.

The most peculiar case found in our work corresponds to a positive static torque for M_p = 0.25 M_J. It is yet to be verified if the torque can remain positive once the planet is allowed to migrate, that is, whether the planet can move upstream in an inward-accreting disk or not.

4.2 Implications for planet formation

In order to fully assess the implications of our findings for the assembly of planetary systems, we would need to derive a migration track for a planet with M_p evolving from a giant-planet core to a fully formed Jupiter (e.g. Bitsch et al. 2015). This is not a straightforward task because we only obtained the static torque, which is not an exact measure of the true migration rate, as pointed out in Sect. 4.1.

Nevertheless, we can at least discuss the migration timescale in a speculative manner, keeping in mind the drawbacks mentioned above. For the cases of M_p = 0.18, 0.25, and 0.5 M_J, which exhibit the most prominent torque reduction, we calculate the migration timescale as (Papaloizou & Larwood 2000) (36)

so that the positive τ_II implies inward migration, and the negative implies outward migration. Without gap irradiation, we obtain τ_II ≃ 45 kyr regardless of M_p. With gap irradiation, we obtain τ_II ≃ 187, − 270, and 157 kyr for M_p = 0.18, 0.25, and 0.5 M_J, respectively. Assuming that the torque would be reduced by an additional factor of 0.2 once the planet is released and a new torque balance is restored (adopted for our parameters from Fig. 14 of Dürmann & Kley 2015), the expected dynamical timescales are then τ_II,dyn ≃ 0.93, − 1.35, and 0.79 Myr. Provided that planets can efficiently form at ≳ 1 au (e.g. Lambrechts & Johansen 2012) and survive Type I migration before the gap opening (e.g. Masset et al. 2006b; Paardekooper & Mellema 2006), we speculate on the basis of τ_II,dyn that Saturn-mass planets with irradiated gaps can easily survive at separations ≳ 1 au while sub-Saturns and half-Jupiters might require special timing for the final stages of their growth with respect to the disk lifetime (τ_disk ≃ 3–10 Myr; Fedele et al. 2010), otherwise they could still substantially migrate inwards. Migration of Jupiter-mass planets can only be slowed down by gap irradiation in other-than-fiducial disks, preferably in those with larger α viscosity.

Our findings can potentially have important implications for the viability of the Grand Tack scenario (Walsh et al. 2011), which critically depends on the migration rate of Saturn. For example, Saturn could never catch up with Jupiter if their migration proceeded as indicated by our Fig. 10. However, our simulations only included one embedded planet and cannot provide any decisive conclusions for a migrating pair of gas giants (two planets would probably cause a more complex temperature perturbation of the irradiated disk).

Finally, we point out that the torque reduction due to gap edge illumination relates in an interesting way to the mechanism of Kanagawa (2019) who proposed that inward migration of giant planets can be stalled or reversed by the dust accumulated at the outer gap edge, which modifies the gap profile via aerodynamic coupling. The mechanism of Kanagawa (2019) becomes stronger with increasing gap width and depth, conversely to what we found for the influence of stellar irradiation. We speculate that both mechanisms can act in a complementary way. If so, we can roughly state that Type II migration is stalled: (i) for Saturn-mass planets due to gap edge illumination; (ii) for Jupiter-mass planets due to dust accumulation; (iii) and for intermediate masses due to a combined contribution of both.

5 Conclusions

Motivated by the results of 2D simulations done by Hallam & Paardekooper (2018), we studied the static torque acting on a gap-openingplanet in a 3D stellar-irradiated passive disk. We investigated how the stellar heating of the exposed outer edge of the planet-induced gap affects the disk-driven torque with the aim of identifying possible consequences for the orbital migration of gas giants.

Our findings can be summarized as follows:

• Gap-opening planets modify the thermal structure of the surrounding disk. Most importantly, they induce a temperature inversion at the outer gap edge. The temperature maps that we obtained are qualitatively similar to the findings of Jang-Condell & Turner (2013);

• The gap depth and width resulting from our 3D radiation hydrodynamics model are similar to predictions based on 2D locally isothermal models (e.g. Kanagawa et al. 2016);

• Vortices are often generated at the outer gap edge and they manifest themselves as fast oscillations of the Type II torque. Our stability analysis suggests that the vortices are possibly excited by the Rossby wave instability of the irradiated gap edge;

• As the outer gap edge becomes puffed up due to more efficient heating, the negative one-sided outer Lindblad torque becomes reduced. However, the effect is less efficient than predicted by Hallam & Paardekooper (2018) because the temperature at the gap edge increases by a factor of ≃ 1.4 in our fully radiative model compared to their estimated increase by a factor of ≃ 4;

• For the viscosity α = 10⁻³ (and other parameters fixed according to Table 1), the total torque is reduced due to gap irradiation in all simulations with M_p < 1 M_J. In summary, the total torque acting on the planet mass M_p = 0.1, 0.18, 0.25, and 0.5 M_J is reduced by a factor of 0.77, 0.22, − 0.17, and 0.3, respectively. The negative reduction factor for M_p = 0.25 M_J (which is close to the mass of Saturn) implies torque reversal, indicating that an outward migration could be possible in this case;

• In the mass range , the torque reduction becomes more prominent with decreasing planet mass. In other words, the reduction becomes strongeras the gap becomes narrower and shallower and the temperature excess recedes towards the planet. For , the reduction becomes weaker with decreasing planet mass because the gap starts to vanish and the temperature excess diminishes accordingly;

• For a Jupiter-mass planet, the torque reduction appears when α > 10⁻³ and becomes stronger with increasing α (again favouring gaps with decreasing width and depth). For the maximum explored value of α = 5 × 10⁻³, the magnitude of the total torque is halved due to gap irradiation.

Our results suggest that the importance of the torque reduction is ultimately determined by an interplay of several competing factors: (i) the behaviour of the Lindblad torque at locations with changing temperature gradients (e.g. Masset 2011); (ii) the gap width (which determines the position and extent of the temperature inversion); (iii) the gapdepth (which determines how much gas is left close to the planet to actually contribute to the torque).

We conclude that the slowdown (or reversal) of inward Type II migration due to gap irradiation is a relevant process, especially for moderately wide and deep gaps. We discussed how it can help to explain the origin of giant planets at moderate separations ≳ au. Since gap irradiation occurs naturally in protoplanetary disks, the effect might play an important role in the assembly of planetary systems.

To provide a final assessment of whether torque reduction is viable, long-term simulations with a mobile planet have to be conducted. The reason is that (i) the static Type II torque is usually larger than the dynamical one (e.g. Scardoni et al. 2020); (ii) the gap-crossing flow and the edge structure (and therefore the influence of irradiation) might be altered if the planet drifts; (iii) the influence of gap edge instabilities on the migrating planet might be different.

Acknowledgements

We wish to thank an anonymous referee whose valuable comments allowed us to significantly improve this paper. The work of OC was supported by The Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations project “e-Infrastructure CZ – LM2018140”.

Appendix A Equation for steady-state radial velocity

Here, we provide details on how to obtain Eq. (18), which is used to calculate the radial velocity v_r after the disk is brought to hydrostatic equilibrium. We start with the general form of the azimuthal component of the momentum equation in a non-rotating reference frame, (A.1)

Let us find a stationary (∂_t = 0) and axially symmetric (∂_θ = 0) form of the equation while also neglecting small terms on the left-hand side (following Takeuchi & Lin 2002). We obtain (A.2)

After expanding the right-hand side and multiplying by rρ, we arrive at (A.3)

which is identical to Eq. (18).

In Eq. (18), the component of the viscous stress tensor, (A.4)

does not depend on v_r once the assumption of ∂_θ = 0 is applied. Similarly, the term dependent on v_ϕ in (A.5)

vanishes. Therefore, Eq. (18) can be used to determine v_r once v_θ, ρ, and ν are known.

Appendix B Comparison to other disk models

The structure of protoplanetary disks is determined by the heating and cooling processes operating within. Our radiative simulations take into account only compressional heating, stellar irradiation, and radiative diffusion and thus it is worthwhile to check how the disk structure would change if the thermal balance was different.

To provide a basic comparison, we calculated two additional equilibrium disks for the fiducial set of parameters (Table 1) and we compared their structure to the case presented in Sect. 3.1. The comparison is shown in Fig. B.1. The first additional disk (red curves) takes into account the viscous heating, which is calculated using the full 3D viscous stress tensor (e.g. Mihalas & Weibel Mihalas 1984) and added to the remaining source terms on the right-hand side of Eq. (5). The second additional disk (blue curves) also adds the viscous heating term but the opacity law is now different; we use the opacity according to Bitsch et al. (2014), which is based on Bell & Lin (1994). The implicationsfor the influence of gap irradiation on Type II migration are discussed in Sect. 4.1.

Fig. B.1 Radial profiles of the midplane temperature T(r) (top) and aspect ratio h(r) (bottom) in an equilibrium disk when passive heating is taken into account (black curve), viscous heating is added (red curve), and the opacity of Bell & Lin (1994) is used instead of that of Flock et al. (2019) (blue curve). One can see that viscous heating increases h(r). In combination with opacity transitions, viscous heating generates bumps in h(r).

References

All Tables

All Figures

Fig. 1 Radial profiles of the accretion mass flux FM (top), aspectratio h (middle), and surface density Σ (bottom) after the hydrostatic relaxation (dashed blue curve) and hydrodynamic relaxation (solid black curve). The latter spanned 6000 Porb. The red dotted curve shows a fit of the h ∝ r2∕7 dependence to the result of the hydrodynamic relaxation. The profiles demonstrate that our equilibrium disk has its global thermodynamics governed by passive heating (stellar irradiation) and its accretion rate is almost uniform.

In the text

Fig. 2 Temperature T0 before planet insertion (top) and the temperature difference T − T0 after gap opening (bottom; taken at t = 3500 Porb) in the vertical plane of the disk. The quantities are azimuthally averaged. The Hill sphere and the location of a Jupiter-mass planet are displayed in the bottom panel. We point out that the scale of the vertical axis is exaggerated. The temperature increase due to the irradiation of the outer gap edge is clearly visible in the bottom panel.

In the text

Fig. 3 Azimuthally averaged profiles of the surface density perturbation Σ∕Σp (top), aspectratio h (middle), and midplane sound speed perturbation cs∕cs,0 (bottom) during the course of the fiducial simulation with irradiation. Solid curves corresponding to various simulation times are distinguished by colour. The vertical dashed line is the planet location, the grey vertical band marks the Hill sphere, and the hatched vertical band shows the extent of the gap width (following the definition of Kanagawa et al. 2016). The horizontal dotted line marks the gap depth according to Kanagawa et al. (2016). Dashed curves in the bottom panel are the least-squares fits of the Gaussian function of Hallam & Paardekooper (2018; magenta) and the skewed Gaussian (turquoise) given by Eq. (27). The peak of h and cs can be directly related to the increased heating of the outer gap edge by stellar irradiation (see Fig. 2).

In the text

Fig. 4 Temporal evolution of the normalized gas-driven static torque γΓ∕Γ0 in the fiducial case. We display the unfiltered measurement with the output sampling of 1∕20 Porb (grey curves), as well as the moving average of the one-sided inner, outer, and total torque (blue, red, and black curve, respectively). The horizontal green dashed line shows the value of the converged torque (averaged over the last 1000 Porb). The final value is γΓ∕Γ0 = −0.035.

In the text

Fig. 5 Radial profile of the specific torque Γ(r) (top) and the cumulative torque Γ(<r) (bottom) for the fiducial set of parameters. The simulation with irradiation (solid black curve) and the reference non-radiative simulation (dotted blue curve) are shown. The grey curves represent 200 samples of Γ(r) recorded over 10 Porb that were used to calculate the black curve as the arithmetic mean and they trace the torque oscillations. The planet location andgap width are shown as in Fig. 3.

In the text

	Fig. 6 Temperature difference in the disk midplane with respect to the unperturbed state (similar to the bottom panel of Fig. 2). The figure is taken from our fiducial simulation at t = 3500 Porb and revealsstreamers (filaments) stretching across the gap and spiral arms, which are offset with respect to the planet-induced perturbations. Some of these structures are marked with arrows.
In the text

	Fig. 7 Surface density perturbation (Σ − Σ0(r))∕Σ0(r) relative to the equilibrium state before planet insertion. The figure is taken at t = 3500 Porb and shows the reference non-radiative simulation (top) and the simulation with stellar irradiation (bottom) for the fiducial set of parameters. The outer gap edge in the bottom panel exhibits additional azimuthal asymmetries.
In the text

	Fig. 8 Midplane vorticity perturbation (ω⊥− ω⊥,0)∕ω⊥,0 relative to the equilibrium state. The figure is taken at t = 3500 Porb and shows the reference non-radiative simulation (top) and the simulation with stellar irradiation (bottom) for the fiducial set of parameters.
In the text

Fig. 9 Azimuthally averaged radial profiles of various hydrodynamic stability criteria in the vicinity of the planet. From top to bottom, we show the Rayleigh criterion (the left-hand side of Eq. (30)), Solberg-Høiland criterion (the left-hand side of Eq. (31)), and the generalized potential vorticity (Eq. (34)). The simulation with gap irradiation (solid black curve) and the reference non-radiative simulation (dotted blue curve) are compared. The grey rectangle shows the extent of the Hill sphere. The inflection point of the black curve at about r ≃ 1.17 rp in the bottom panel suggests that the disk is Rossby-unstable.

In the text

Fig. 10 Dependence of the normalized torque γΓ∕Γ0 on the planet mass Mp (top) and α viscosity (bottom). Additionally, the top horizontal axis of each panel shows the position of the outer gap edge redge as resulting from Eq. (24). Circles and solid black lines correspond to simulations with gap irradiation; diamonds and dashed grey lines correspond to reference non-radiative simulations. Open symbols correspond to simulations without torque oscillations (i.e. Rossby-stable) while filled symbols mark simulations with gap edge instabilities. Data points are labelled with the corresponding torque reduction factor firr (Eq. (35)). Clearly, gap edge irradiation reduces the magnitude of the torque (or even causes torque reversal) for Mp < 1 MJ and α > 10−3.

In the text

Fig. 11 Radial profiles of the surface density perturbation Σ∕Σp (top), temperature T (middle), and specific torque Γ(r) (bottom). Simulations with different planet masses Mp = 0.25 (red), 0.5 (blue), and 1 MJ (black) are shown. Simulations with irradiation and reference non-radiative simulations are distinguished using solid and dotted lines,respectively. The top two panels are azimuthally averaged.

In the text

Fig. 12 Specific torque Γ(r) as a function of the radial separation from the planet for various planetary masses (as in Fig. 11). The predominantly positive branch of the curves corresponds to the inner one-sided torque (arising from the interval 0.65 rp < r < 1 rp) and the negative branch is the outer one-sided torque (arising from the interval 1 rp < r< 1.35 rp). The aim of the figure is to assess the inner-outer asymmetry of the one-sided torques and Γ(r) is therefore artificially normalized to the maximum value of the inner peak. The figure reveals how gap irradiation reduces the outer Lindblad torque with respect to the inner one.

In the text

	Fig. 13 Specific cumulative torque Γ(< r) as a function of the radial distance for various planetary masses (as in Fig. 11). Simulations with gap irradiation (top) and reference non-radiative simulations (bottom) are compared. The figure demonstrates how gap irradiation modifies the drop of the cumulative torque across the outer half of the gap.
In the text

	Fig. 14 As in Fig. 12 but for simulations with different values of the α viscosity 10−3 (black), 3 × 10−3 (blue), and 5 × 10−3 (red).
In the text

Fig. B.1 Radial profiles of the midplane temperature T(r) (top) and aspect ratio h(r) (bottom) in an equilibrium disk when passive heating is taken into account (black curve), viscous heating is added (red curve), and the opacity of Bell & Lin (1994) is used instead of that of Flock et al. (2019) (blue curve). One can see that viscous heating increases h(r). In combination with opacity transitions, viscous heating generates bumps in h(r).

In the text