Issue 
A&A
Volume 642, October 2020



Article Number  A219  
Number of page(s)  17  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202038988  
Published online  23 October 2020 
Migration of gapopening planets in 3D stellarirradiated accretion disks
^{1}
Astronomical Institute of Charles University,
V Holešovičkách 747/2,
180 00
Prague 8,
Czech Republic
email: chrenko@sirrah.troja.mff.cuni.cz
^{2}
Department of Space Studies, Southwest Research Institute,
1050 Walnut Street, Suite 300,
Boulder,
CO
80302,
USA
Received:
21
July
2020
Accepted:
15
September
2020
Context. The origin of giant planets at moderate separations ≃1–10 au is still not fully understood because numerical studies of Type II migration in protoplanetary disks often predict a decay of the semimajor axis that is too fast. According to recent 2D simulations, inward migration of a gapopening planet can be slowed down or even reversed if the outer gap edge becomes heated by irradiation from the central star, and puffed up.
Aims. Here, we study how stellar irradiation reduces the diskdriven torque and affects migration in more realistic 3D disks.
Methods. Using 3D hydrodynamic simulations with radiation transfer, we investigated the static torque acting on a single gapopening planet embedded in a passively heated accretion disk.
Results. Our simulations confirm that a temperature inversion is established at the irradiated outer gap edge and the local increase of the scale height reduces the magnitude of the negative outer Lindblad torque. However, the temperature excess is smaller than assumed in 2D simulations and the torque reduction only becomes prominent for specific parameters. For the viscosity α = 10^{−3}, the total torque is reduced for planetary masses ranging from 0.1 to 0.7 Jupiter mass, with the strongest reduction being by a factor of − 0.17 (implying outward migration) for a Saturnmass planet. For a Jupitermass planet, the torque reduction becomes stronger with increasing α (the torque is halved when α = 5 × 10^{−3}).
Conclusions. We conclude that planets that open moderately wide and deep gaps are subject to the largest torque modifications and their Type II migration can be stalled due to gap edge illumination. We then argue that the torque reduction can help to stabilize the orbits of giant planets forming at ≳ 1 au.
Key words: hydrodynamics / planets and satellites: formation / planetdisk interactions / protoplanetary disks
© ESO 2020
1 Introduction
Giant planets form in protoplanetary disks and their early evolution is driven by gravitational planetdisk 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 razorthin accretion disks (Shakura & Sunyaev 1973; LyndenBell & 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 longterm 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^{3} 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 onesided Lindblad torque scales as Γ_{OS} ~ Ch^{−3} (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 gapopening 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 nonisothermal 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ítezLlambay & 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 fluxlimited 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 stellarirradiation term Q_{irr} while the heating by viscous dissipation is not considered^{1}. 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 selfshadowing 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 selfshadowing can only occur once the planet is inserted and starts to perturb the gas distribution.
Stellar irradiation is implemented following DobbsDixon 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 τ_{0} 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 crosssection of the cell, and V_{cell} is its volume.
Regarding the disk opacity, we assume it is dominated by submicron 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^{2} g^{−1} where κ_{R} is the Rosseland mean opacity (which enters the model through the FLD). The opacity to stellar irradiation is κ_{⋆} = 1300 cm^{2} g^{−1}. To evaluate theopacity of a gas–dust mixture, each opacity value is scaled by the dusttogas ratio, which we choose as Z = 0.001. The somewhat lower value of Z reflects the fact that submicron grains tend to get depleted as the mass spectrum of the coagulationfragmentation 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). Dustfree regions would have a reduced opacity κ_{gas} ≃ 10^{−5} cm^{2} g^{−1} (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 temperaturedependent 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 righthand side.
The gravitational potential generated by the star and the planet is (10)
where d is the cellplanet 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 selfgravity 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 diskdriven torque. The cutoff function is (Crida et al. 2008) (11)
and we use p = 0.6 (Robert et al. 2018).
The twobody starplanet 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; .
Summary of the fiducial parameters.
2.2 Parameters
The parameters that we use in our fiducial simulation are listed in Table 1. We focus on solartype 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 Jupitermass 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^{−8} M_{⊙} yr^{−1} is close to the mean value observed around solarmass protostars in Lupus and Chamaeleon I (Manara et al. 2017; Mulders et al. 2017). The fiducial viscosity α = 10^{−3} is applicableto disks in which the angular momentum transport is facilitated by the hydrodynamic turbulence or disk winds (α ≃ 10^{−3}–10^{−4}; Nelson et al. 2013; Klahr & Hubbard 2014; Béthune et al. 2017) rather than magnetorotational instability (α ≃ 10^{−1}–10^{−2}; 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 planetdisk 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 nonisothermal 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 fourstage 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 planetdisk 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 razorthin 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 τ_{0}(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ítezLlambay 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 overrelaxation method with the relative precision ɛ_{sor} = 10^{−8}. 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 densityweighted vertically averaged viscosity (16)
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^{−5} 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 powerlaw 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_{V}T in each ghost cell.
Additionally, the boundary conditions are supplemented with the wavekilling zones of de ValBorro 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 τ_{0} (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 nonaccreting (K = 0), unless stated otherwise.
2.3.4 Main stage
After the planet insertion, the simulation is continued until the measured diskdriven 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^{6}); (ii) fast wave propagation in lowdensity 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 wavekilling procedure and by introducing a volume density floor ρ_{floor} = 10^{−22} g cm^{−3}.
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 MultiProcessing (OpenMP). A single simulation was usually spawned over ≃ 550 CPU cores and consumed ~ 10^{5} CPU hours.
2.4 Reference nonradiative 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^{−3} 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 F_{M} (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 P_{orb}. The red dotted curve shows a fit of the h ∝ r^{2∕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), planetinduced 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^{−8} M_{⊙} yr^{−1}. 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 leastsquares fit of h ∝ r^{2∕7} in accordance with Chiang & Goldreich (1997). The Σ(r) profile exhibits Σ_{p} = 385 g cm^{−2} 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 T_{0} before planet insertion (top) and the temperature difference T − T_{0} after gap opening (bottom; taken at t = 3500 P_{orb}) in the vertical plane of the disk. The quantities are azimuthally averaged. The Hill sphere and the location of a Jupitermass 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 largescale 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 JangCondell & 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)
By plugging in the values from our simulation, we find the gap depth to be Σ_{min}∕Σ_{0} ≃ 1.2 × 10^{−3}, 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^{2} 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 puffedup 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 leastsquares 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 leastsquares 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 c_{s}∕c_{s,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 leastsquares fits of the Gaussian function of Hallam & Paardekooper (2018; magenta) and the skewed Gaussian (turquoise) given by Eq. (27). The peak of h and c_{s} 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 gasdriven static torque γΓ∕Γ_{0} in the fiducial case. We display the unfiltered measurement with the output sampling of 1∕20 P_{orb} (grey curves), as well as the moving average of the onesided 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 P_{orb}). 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 γΓ∕Γ_{0} where (e.g. Paardekooper et al. 2010) (28)
and all quantities correspond to the state before planet insertion. The quantity Γ_{0} is defined to reflect the basic dependencies (e.g. ; ~ h^{−2}) 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 γΓ∕Γ_{0}, the factor γ accounts forthe difference in the sound speed between isothermal and nonisothermal models (Baruteau & Masset 2008).
From Fig. 4, one can see that the total torque undergoes fast oscillations related to the outer onesided 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} = −0.035.
The question now arises – is the measured torque reduced due to gap edge irradiation? To answer the question, we performed a reference nonradiative 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^{3}.
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 nonradiative simulation (dotted blue curve) are shown. The grey curves represent 200 samples of Γ(r) recorded over 10 P_{orb} 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 P_{orb} and revealsstreamers (filaments) stretching across the gap and spiral arms, which are offset with respect to the planetinduced 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 coorbital region or at gap edges (Koller et al. 2003; Li et al. 2005; de ValBorro 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 wavelike 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 nonzero 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 coorbital 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 P_{orb} and shows the reference nonradiative 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 P_{orb} and shows the reference nonradiative 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 planetinduced 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^{2}Ω increases with radius (Chandrasekhar 1961). For the Rayleigh stability criterion, we use (30)
where the brackets denote the vertical densityweighted average (as in Eq. (16)).
The buoyant instability appears in nonbarotropic 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 planetdriven shocks (e.g. Ou et al. 2007; Richert et al. 2015). To assess the buoyant stability, we adopt a form of the SolbergHø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 BruntVä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 planetinduced gaps. It can become excited when the generalized potential vorticity (Lin 2013) (34)
develops a local inflection point. The quantity ξ essentially represents an entropymodified 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 lefthand side of Eq. (30)), SolbergHøiland criterion (the lefthand side of Eq. (31)), and the generalized potential vorticity (Eq. (34)). The simulation with gap irradiation (solid black curve) and the reference nonradiative 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 r_{p} in the bottom panel suggests that the disk is Rossbyunstable. 
3.2 Dependence on parameters
In this section, we perform a coarse parametric study to test the behaviour of the Type II torque in stellarirradiated 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^{−9} M_{⊙} yr^{−1}, 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^{−4}, 5 × 10^{−3}, 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^{−3}. 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 nonradiative 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^{−3}, 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^{−3} (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^{−3}. 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 γΓ∕Γ_{0}(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 γΓ∕Γ_{0}(α) 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^{−9} M_{⊙} yr^{−1}, the planet is considerably more massive than the disk (M_{D} = 0.12 M_{J}) and consequently, the planetinduced perturbations of the outer disk are so strong that neither the disk structure nor the torque can converge before 5000 P_{orb}. For Ṁ = 10^{−7} M_{⊙} yr^{−1}, 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^{−9} and 5 × 10^{−8} M_{⊙} yr^{−1} proceeded normally and we obtained the normalized value of the torque as γΓ∕Γ_{0} = −0.032 and − 0.037, respectively.Both these values are very similar to the fiducial simulation with γΓ∕Γ_{0} = −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.
Summary of our parametric study.
Fig. 10
Dependence of the normalized torque γΓ∕Γ_{0} on the planet mass M_{p} (top) and α viscosity (bottom). Additionally, the top horizontal axis of each panel shows the position of the outer gap edge r_{edge} as resulting from Eq. (24). Circles and solid black lines correspond to simulations with gap irradiation; diamonds and dashed grey lines correspond to reference nonradiative simulations. Open symbols correspond to simulations without torque oscillations (i.e. Rossbystable) while filled symbols mark simulations with gap edge instabilities. Data points are labelled with the corresponding torque reduction factor f_{irr} (Eq. (35)). Clearly, gap edge irradiation reduces the magnitude of the torque (or even causes torque reversal) for M_{p} < 1 M_{J} 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 M_{p} = 0.25 (red), 0.5 (blue), and 1 M_{J} (black) are shown. Simulations with irradiation and reference nonradiative 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 torquemass 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 onesided torques.
As a confirmation, we compare the onesided 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 nonradiative simulations, the innerouter asymmetry is always such that the outer negative onesided torque dominates. However, in simulations with gap irradiation and M_{p} = 0.25 and 0.5 M_{J}, the extremal values of the onesided torques are almost identical. For M_{p} = 0.5 M_{J}, the outer onesided 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 onesided torques in the depicted interval of Δr almost cancel out. Interestingly, we notice that the innerouter 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 onesided torque (arising from the interval 0.65 r_{p} < r < 1 r_{p}) and the negative branch is the outer onesided torque (arising from the interval 1 r_{p} < r< 1.35 r_{p}). The aim of the figure is to assess the innerouter asymmetry of the onesided 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 nonradiative 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 innerouter asymmetry of the onesided torques for the fiducial case (α = 10^{−3}) and cases α = 3 × 10^{−3} and 5 × 10^{−3}. 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^{−3}, we identify that (i) the innerouter 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 innerouter asymmetry of the onesided 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 semimajor 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}∕Σ_{0} (Eq. (25)) decreases with increasing r. Therefore, our results cannot be directly scaled to an arbitrary planetary semimajor axis and some trends can be expected. For example, the torque acting on a Jupitermass 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 Σ_{0}) and the temperature excess would be shifted as well. The Jupitermass 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 selfshadowed 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 innerouter asymmetry of the onesided 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, gapcrossing 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 inwardaccreting 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 giantplanet 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 Saturnmass planets with irradiated gaps can easily survive at separations ≳ 1 au while subSaturns and halfJupiters 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 Jupitermass planets can only be slowed down by gap irradiation in otherthanfiducial 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 Saturnmass planets due to gap edge illumination; (ii) for Jupitermass 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 gapopeningplanet in a 3D stellarirradiated passive disk. We investigated how the stellar heating of the exposed outer edge of the planetinduced gap affects the diskdriven torque with the aim of identifying possible consequences for the orbital migration of gas giants.
Our findings can be summarized as follows:

Gapopening 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 JangCondell & 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 onesided 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^{−3} (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 Jupitermass planet, the torque reduction appears when α > 10^{−3} and becomes stronger with increasing α (again favouring gaps with decreasing width and depth). For the maximum explored value of α = 5 × 10^{−3}, 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, longterm 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 gapcrossing 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 “eInfrastructure CZ – LM2018140”.
Appendix A Equation for steadystate 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 nonrotating 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 lefthand side (following Takeuchi & Lin 2002). We obtain (A.2)
After expanding the righthand 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 righthand 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
 Baraffe, I., Chabrier, G., Allard, F., & Hauschildt, P. H. 1998, A&A, 337, 403 [Google Scholar]
 Baruteau, C., & Masset, F. 2008, ApJ, 672, 1054 [NASA ADS] [CrossRef] [Google Scholar]
 Bell, K. R., & Lin, D. N. C. 1994, ApJ, 427, 987 [NASA ADS] [CrossRef] [Google Scholar]
 BenítezLlambay, P., & Masset, F. S. 2016, ApJS, 223, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., & Kley, W. 2010, A&A, 523, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Crida, A., Morbidelli, A., Kley, W., & DobbsDixon, I. 2013, A&A, 549, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Morbidelli, A., Lega, E., & Crida, A. 2014, A&A, 564, A135 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Lambrechts, M., & Johansen, A. 2015, A&A, 582, A112 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bitsch, B., Izidoro, A., Johansen, A., et al. 2019, A&A, 623, A88 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cassan, A., Kubas, D., Beaulieu, J. P., et al. 2012, Nature, 481, 167 [NASA ADS] [CrossRef] [Google Scholar]
 Chandrasekhar, S. 1961, Hydrodynamic and Hydromagnetic Stability (Oxford: Clarendon Press) [Google Scholar]
 Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 [NASA ADS] [CrossRef] [Google Scholar]
 Chrenko, O., & Lambrechts, M. 2019, A&A, 626, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Coleman, G. A. L., & Nelson, R. P. 2014, MNRAS, 445, 479 [NASA ADS] [CrossRef] [Google Scholar]
 Coleman, G. A. L., & Nelson, R. P. 2016, MNRAS, 460, 2779 [NASA ADS] [CrossRef] [Google Scholar]
 Crida, A., & Bitsch, B. 2017, Icarus, 285, 145 [NASA ADS] [CrossRef] [Google Scholar]
 Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587 [Google Scholar]
 Crida, A., Sándor, Z., & Kley, W. 2008, A&A, 483, 325 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Crida, A., Bitsch, B., & Raibaldi, A. 2016, in SF2A2016: Proceedings of the Annual Meeting of the French Society of Astronomy and Astrophysics, eds. C. Reylé, J. Richard, L. Cambrésy, et al., 473 [Google Scholar]
 de ValBorro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529 [NASA ADS] [CrossRef] [Google Scholar]
 de ValBorro, M., Artymowicz, P., D’Angelo, G., & Peplinski, A. 2007, A&A, 471, 1043 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 DobbsDixon, I., Cumming, A., & Lin, D. N. C. 2010, ApJ, 710, 1395 [NASA ADS] [CrossRef] [Google Scholar]
 Duffell, P. C., Haiman, Z., MacFadyen, A. I., D’Orazio, D. J., & Farris, B. D. 2014, ApJ, 792, L10 [NASA ADS] [CrossRef] [Google Scholar]
 Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957 [Google Scholar]
 Dürmann, C., & Kley, W. 2015, A&A, 574, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dzyurkevich, N., Turner, N. J., Henning, T., & Kley, W. 2013, ApJ, 765, 114 [NASA ADS] [CrossRef] [Google Scholar]
 Fedele, D., van den Ancker, M. E., Henning, T., Jayawardhana, R., & Oliveira, J. M. 2010, A&A, 510, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Flock, M., Nelson, R. P., Turner, N. J., et al. 2017, ApJ, 850, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Flock, M., Turner, N. J., Mulders, G. D., et al. 2019, A&A, 630, A147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics: Third Edition (Cambridge, UK: Cambridge University Press) [Google Scholar]
 Fressin, F., Torres, G., Charbonneau, D., et al. 2013, ApJ, 766, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Fromang, S., & Nelson, R. P. 2006, A&A, 457, 343 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Fromang, S., Lyra, W., & Masset, F. 2011, A&A, 534, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fu, W., Li, H., Lubow, S., & Li, S. 2014, ApJ, 788, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Fung, J., & Chiang, E. 2016, ApJ, 832, 105 [NASA ADS] [CrossRef] [Google Scholar]
 Hallam, P. D., & Paardekooper, S. J. 2018, MNRAS, 481, 1667 [NASA ADS] [CrossRef] [Google Scholar]
 Hasegawa, Y., & Ida, S. 2013, ApJ, 774, 146 [NASA ADS] [CrossRef] [Google Scholar]
 Isella, A., & Natta, A. 2005, A&A, 438, 899 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jacquet, E. 2013, A&A, 551, A75 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 JangCondell, H., & Turner, N. J. 2013, ApJ, 772, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Johansen, A., Ida, S., & Brasser, R. 2019, A&A, 622, A202 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kanagawa, K. D. 2019, ApJ, 879, L19 [CrossRef] [Google Scholar]
 Kanagawa, K. D., Tanaka, H., Muto, T., Tanigawa, T., & Takeuchi, T. 2015, MNRAS, 448, 994 [NASA ADS] [CrossRef] [Google Scholar]
 Kanagawa, K. D., Muto, T., Tanaka, H., et al. 2016, PASJ, 68, 43 [NASA ADS] [Google Scholar]
 Kanagawa, K. D., Tanaka, H., & Szuszkiewicz, E. 2018, ApJ, 861, 140 [NASA ADS] [CrossRef] [Google Scholar]
 Klahr, H. 2004, ApJ, 606, 1070 [Google Scholar]
 Klahr, H., & Hubbard, A. 2014, ApJ, 788, 21 [Google Scholar]
 Klahr, H., & Kley, W. 2006, A&A, 445, 747 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Klahr, H. H., & Bodenheimer, P. 2003, ApJ, 582, 869 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W. 1989, A&A, 208, 98 [NASA ADS] [Google Scholar]
 Kley, W. 1999, MNRAS, 303, 696 [NASA ADS] [CrossRef] [Google Scholar]
 Kley, W., D’Angelo, G., & Henning, T. 2001, ApJ, 547, 457 [NASA ADS] [CrossRef] [Google Scholar]
 Kolb, S. M., Stute, M., Kley, W., & Mignone, A. 2013, A&A, 559, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Koller, J., Li, H., & Lin, D. N. C. 2003, ApJ, 596, L91 [NASA ADS] [CrossRef] [Google Scholar]
 Lambrechts, M., & Johansen, A. 2012, A&A, 544, A32 [Google Scholar]
 Les, R., & Lin, M.K. 2015, MNRAS, 450, 1503 [NASA ADS] [CrossRef] [Google Scholar]
 Lesur, G., & Papaloizou, J. C. B. 2010, A&A, 513, A60 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
 Li, H., Finn, J. M., Lovelace, R. V. E., & Colgate, S. A. 2000, ApJ, 533, 1023 [NASA ADS] [CrossRef] [Google Scholar]
 Li, H., Li, S., Koller, J., et al. 2005, ApJ, 624, 1003 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, M.K. 2013, ApJ, 765, 84 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. 1986a, ApJ, 307, 395 [Google Scholar]
 Lin, D. N. C., & Papaloizou, J. 1986b, ApJ, 309, 846 [NASA ADS] [CrossRef] [Google Scholar]
 Lin, M.K., & Papaloizou, J. C. B. 2010, MNRAS, 405, 1473 [NASA ADS] [Google Scholar]
 Lovelace, R. V. E., Li, H., Colgate, S. A., & Nelson, A. F. 1999, ApJ, 513, 805 [NASA ADS] [CrossRef] [Google Scholar]
 Lubow, S. H., & D’Angelo, G. 2006, ApJ, 641, 526 [NASA ADS] [CrossRef] [Google Scholar]
 LyndenBell, D., & Pringle, J. E. 1974, MNRAS, 168, 603 [NASA ADS] [CrossRef] [Google Scholar]
 Manara, C. F., Testi, L., Herczeg, G. J., et al. 2017, A&A, 604, A127 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Masset, F. S. 2011, Celes. Mech. Dyn. Astron., 111, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Masset, F. S., & BenítezLlambay, P. 2016, ApJ, 817, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Masset, F. S., D’Angelo, G., & Kley, W. 2006a, ApJ, 652, 730 [NASA ADS] [CrossRef] [Google Scholar]
 Masset, F. S., Morbidelli, A., Crida, A., & Ferreira, J. 2006b, ApJ, 642, 478 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumura, S., & Pudritz, R. E. 2005, ApJ, 618, L137 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., Marmier, M., Lovis, C., et al. 2011, ArXiv eprints, [arXiv:1109.2497] [Google Scholar]
 Meheut, H., Casse, F., Varniere, P., & Tagger, M. 2010, A&A, 516, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Mihalas, D., & Weibel Mihalas, B. 1984, Foundations of Radiation Hydrodynamics (New York: Oxford University Press) [Google Scholar]
 Mulders, G. D., Pascucci, I., Manara, C. F., et al. 2017, ApJ, 847, 31 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, R. P., Papaloizou, J. C. B., Masset, F., & Kley, W. 2000, MNRAS, 318, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Nelson, R. P., Gressel, O., & Umurhan, O. M. 2013, MNRAS, 435, 2610 [NASA ADS] [CrossRef] [Google Scholar]
 Ou, S., Ji, J., Liu, L., & Peng, X. 2007, ApJ, 667, 1220 [NASA ADS] [CrossRef] [Google Scholar]
 Paardekooper,S.J., & Mellema, G. 2006, A&A, 459, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Paardekooper,S.J., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401, 1950 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., & Larwood, J. D. 2000, MNRAS, 315, 823 [NASA ADS] [CrossRef] [Google Scholar]
 Papaloizou, J. C. B., Nelson, R. P., Kley, W., Masset, F. S., & Artymowicz, P. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil, 655 [Google Scholar]
 Petersen, M. R., Julien, K., & Stewart, G. R. 2007, ApJ, 658, 1236 [NASA ADS] [CrossRef] [Google Scholar]
 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Rein, H., & Liu, S.F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424 [NASA ADS] [CrossRef] [Google Scholar]
 Richert, A. J. W., Lyra, W., Boley, A., Mac Low, M.M., & Turner, N. 2015, ApJ, 804, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Robert, C. M. T., Crida, A., Lega, E., Méheut, H., & Morbidelli, A. 2018, A&A, 617, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rüdiger, G., Arlt, R., & Shalybkov, D. 2002, A&A, 391, 781 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Santerne, A., Moutou, C., Tsantaki, M., et al. 2016, A&A, 587, A64 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Scardoni, C. E., Rosotti, G. P., Lodato, G., & Clarke, C. J. 2020, MNRAS, 492, 1318 [CrossRef] [Google Scholar]
 Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 500, 33 [NASA ADS] [Google Scholar]
 Stone, J. M., & Norman, M. L. 1992, ApJS, 80, 753 [Google Scholar]
 Takeuchi, T., & Lin, D. N. C. 2002, ApJ, 581, 1344 [NASA ADS] [CrossRef] [Google Scholar]
 Terquem, C. E. J. M. L. J. 2008, ApJ, 689, 532 [NASA ADS] [CrossRef] [Google Scholar]
 Walsh, K. J., Morbidelli, A., Raymond, S. N., O’Brien, D. P., & Mandell, A. M. 2011, Nature, 475, 206 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Ward, W. R. 1986, Icarus, 67, 164 [Google Scholar]
 Ward, W. R. 1997, Icarus, 126, 261 [Google Scholar]
 White, R. J., Greene, T. P., Doppmann, G. W., Covey, K. R., & Hillenbrand, L. A. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil, 117 [Google Scholar]
 Zhu, Z., Dong, R., Stone, J. M., & Rafikov, R. R. 2015, ApJ, 813, 88 [NASA ADS] [CrossRef] [Google Scholar]
 Ziampras, A., Kley, W., & Dullemond, C. P. 2020, A&A, 637, A50 [Google Scholar]
We point out, however, that the artificial viscosity term of FARGO3D, which is used for spreading shocks (Stone & Norman 1992), is included and the resulting heating term is accounted for.
We point out that the definition of the gap width operates with the ratio Σ∕Σ_{0}(r) while Fig. 3 displays Σ∕Σ_{p} (Σ_{0}(r) is the unperturbed surface density at a given radius; Σ_{p} is the unperturbed surface density at the planet location). We derived Δ_{gap} using the correct quantity.
One might even argue that gap irradiation increases the magnitude of the negative total torque in this case but we believe that the small difference, ΔγΓ∕Γ_{0} = 0.004, might also be attributed to the differences in the physical treatment of the energy used in our model with irradiation and reference nonradiative model.
All Tables
All Figures
Fig. 1
Radial profiles of the accretion mass flux F_{M} (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 P_{orb}. The red dotted curve shows a fit of the h ∝ r^{2∕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 T_{0} before planet insertion (top) and the temperature difference T − T_{0} after gap opening (bottom; taken at t = 3500 P_{orb}) in the vertical plane of the disk. The quantities are azimuthally averaged. The Hill sphere and the location of a Jupitermass 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 c_{s}∕c_{s,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 leastsquares fits of the Gaussian function of Hallam & Paardekooper (2018; magenta) and the skewed Gaussian (turquoise) given by Eq. (27). The peak of h and c_{s} 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 gasdriven static torque γΓ∕Γ_{0} in the fiducial case. We display the unfiltered measurement with the output sampling of 1∕20 P_{orb} (grey curves), as well as the moving average of the onesided 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 P_{orb}). 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 nonradiative simulation (dotted blue curve) are shown. The grey curves represent 200 samples of Γ(r) recorded over 10 P_{orb} 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 P_{orb} and revealsstreamers (filaments) stretching across the gap and spiral arms, which are offset with respect to the planetinduced 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 P_{orb} and shows the reference nonradiative 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 P_{orb} and shows the reference nonradiative 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 lefthand side of Eq. (30)), SolbergHøiland criterion (the lefthand side of Eq. (31)), and the generalized potential vorticity (Eq. (34)). The simulation with gap irradiation (solid black curve) and the reference nonradiative 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 r_{p} in the bottom panel suggests that the disk is Rossbyunstable. 

In the text 
Fig. 10
Dependence of the normalized torque γΓ∕Γ_{0} on the planet mass M_{p} (top) and α viscosity (bottom). Additionally, the top horizontal axis of each panel shows the position of the outer gap edge r_{edge} as resulting from Eq. (24). Circles and solid black lines correspond to simulations with gap irradiation; diamonds and dashed grey lines correspond to reference nonradiative simulations. Open symbols correspond to simulations without torque oscillations (i.e. Rossbystable) while filled symbols mark simulations with gap edge instabilities. Data points are labelled with the corresponding torque reduction factor f_{irr} (Eq. (35)). Clearly, gap edge irradiation reduces the magnitude of the torque (or even causes torque reversal) for M_{p} < 1 M_{J} 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 M_{p} = 0.25 (red), 0.5 (blue), and 1 M_{J} (black) are shown. Simulations with irradiation and reference nonradiative 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 onesided torque (arising from the interval 0.65 r_{p} < r < 1 r_{p}) and the negative branch is the outer onesided torque (arising from the interval 1 r_{p} < r< 1.35 r_{p}). The aim of the figure is to assess the innerouter asymmetry of the onesided 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 nonradiative 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.