Oscillatory migration of accreting protoplanets driven by a 3D   distortion of the gas flow

Ond\v{r}ej Chrenko; Michiel Lambrechts

arXiv:1904.12497·astro-ph.EP·June 26, 2019

Oscillatory migration of accreting protoplanets driven by a 3D distortion of the gas flow

Ond\v{r}ej Chrenko, Michiel Lambrechts

PDF

TL;DR

This study uses 3D radiation hydrodynamic simulations to show that temperature-dependent opacities cause oscillatory migration of accreting protoplanets due to flow instabilities and asymmetries in the gas flow.

Contribution

It reveals how non-uniform opacities induce flow instabilities that lead to oscillatory torques and migration patterns in protoplanets, extending previous steady-state models.

Findings

01

Heating torque causes flow perturbations and asymmetries.

02

Flow instabilities lead to oscillatory protoplanet migration.

03

Oscillations result in alternating inward and outward movement.

Abstract

Context. The dynamics of a low-mass protoplanet accreting solids is influenced by the heating torque, which was found to suppress inward migration in protoplanetary disks with constant opacities. Aims. We investigate the differences of the heating torque between disks with constant and temperature-dependent opacities. Methods. Interactions of a super-Earth-sized protoplanet with the gas disk are explored using 3D radiation hydrodynamic simulations. Results. Accretion heating of the protoplanet creates a hot underdense region in the surrounding gas, leading to misalignment of the local density and pressure gradients. As a result, the 3D gas flow is perturbed and some of the streamlines form a retrograde spiral rising above the protoplanet. In the constant-opacity disk, the perturbed flow reaches a steady state and the underdense gas responsible for the heating torque remains…

Equations60

\frac{\partial ρ}{\partial t} + (\@vec v \cdot \nabla) ρ = - ρ \nabla \cdot \@vec v,

\frac{\partial ρ}{\partial t} + (\@vec v \cdot \nabla) ρ = - ρ \nabla \cdot \@vec v,

\frac{\partial \@vec v}{\partial t} + (\@vec v \cdot \nabla) \@vec v = - \frac{\nabla P}{ρ} - \nablaΦ + \frac{\nabla \cdot \@tens T}{ρ} - [2 \@vec Ω \times \@vec v + \@vec Ω \times (\@vec Ω \times \@vec r)],

\frac{\partial \@vec v}{\partial t} + (\@vec v \cdot \nabla) \@vec v = - \frac{\nabla P}{ρ} - \nablaΦ + \frac{\nabla \cdot \@tens T}{ρ} - [2 \@vec Ω \times \@vec v + \@vec Ω \times (\@vec Ω \times \@vec r)],

\frac{\partial E _{R}}{\partial t} + \nabla \cdot \@vec F = ρ κ_{P} [4 σ T^{4} - c E_{R}],

\frac{\partial E _{R}}{\partial t} + \nabla \cdot \@vec F = ρ κ_{P} [4 σ T^{4} - c E_{R}],

\frac{\partial ϵ}{\partial t} + (\@vec v \cdot \nabla) ϵ = - P \nabla \cdot \@vec v - ρ κ_{P} [4 σ T^{4} - c E_{R}] + Q_{visc} + Q_{art} + Q_{acc},

\frac{\partial ϵ}{\partial t} + (\@vec v \cdot \nabla) ϵ = - P \nabla \cdot \@vec v - ρ κ_{P} [4 σ T^{4} - c E_{R}] + Q_{visc} + Q_{art} + Q_{acc},

P = (γ - 1) ϵ = (γ - 1) ρ c_{V} T,

P = (γ - 1) ϵ = (γ - 1) ρ c_{V} T,

\@vec F = - λ_{lim} \frac{c}{ρ κ _{R}} \nabla E_{R},

\@vec F = - λ_{lim} \frac{c}{ρ κ _{R}} \nabla E_{R},

L = \frac{G M _{p} M ˙ _{p}}{R _{p}} = \frac{G M _{p}^{2}}{R _{p} τ},

L = \frac{G M _{p} M ˙ _{p}}{R _{p}} = \frac{G M _{p}^{2}}{R _{p} τ},

Q_{acc} = \frac{L}{V},

Q_{acc} = \frac{L}{V},

Φ = - \frac{G M _{⋆}}{r} - \frac{G M _{p}}{d} f_{sm},

Φ = - \frac{G M _{⋆}}{r} - \frac{G M _{p}}{d} f_{sm},

f_{\mathrm{sm}}=\begin{dcases*}1\,,&$\left(d>r_{\mathrm{sm}}\right)$ ,\\ \left[\left(\frac{d}{r_{\mathrm{sm}}}\right)^{4}-2\left(\frac{d}{r_{\mathrm{sm}}}\right)^{3}+2\frac{d}{r_{\mathrm{sm}}}\right]\,,&$\left(d\leq r_{\mathrm{sm}}\right)$\end{dcases*}

f_{\mathrm{sm}}=\begin{dcases*}1\,,&$\left(d>r_{\mathrm{sm}}\right)$ ,\\ \left[\left(\frac{d}{r_{\mathrm{sm}}}\right)^{4}-2\left(\frac{d}{r_{\mathrm{sm}}}\right)^{3}+2\frac{d}{r_{\mathrm{sm}}}\right]\,,&$\left(d\leq r_{\mathrm{sm}}\right)$\end{dcases*}

Γ = disk \int (\@vec r_{p} \times \@vec F_{g})_{⊥} d V,

Γ = disk \int (\@vec r_{p} \times \@vec F_{g})_{⊥} d V,

λ_{c} = \frac{χ}{q Ω _{p} γ},

λ_{c} = \frac{χ}{q Ω _{p} γ},

\@vec ω_{a} \equiv \@vec ω_{r} + 2 \@vec Ω \equiv \nabla \times \@vec v + 2 \@vec Ω,

\@vec ω_{a} \equiv \@vec ω_{r} + 2 \@vec Ω \equiv \nabla \times \@vec v + 2 \@vec Ω,

\frac{D \@vec ω _{r}}{D t} = (\@vec ω_{a} \cdot \nabla) \@vec v - \@vec ω_{a} (\nabla \cdot \@vec v) + \frac{\nabla ρ \times \nabla P}{ρ ^{2}},

\frac{D \@vec ω _{r}}{D t} = (\@vec ω_{a} \cdot \nabla) \@vec v - \@vec ω_{a} (\nabla \cdot \@vec v) + \frac{\nabla ρ \times \nabla P}{ρ ^{2}},

∣ \nabla_{rad, ϕ} ∣ > ∣ \nabla_{ad, ϕ} ∣ \Leftrightarrow \frac{∣ \nabla _{rad, ϕ} ∣}{∣ \nabla _{ad, ϕ} ∣} - 1 > 0,

∣ \nabla_{rad, ϕ} ∣ > ∣ \nabla_{ad, ϕ} ∣ \Leftrightarrow \frac{∣ \nabla _{rad, ϕ} ∣}{∣ \nabla _{ad, ϕ} ∣} - 1 > 0,

(\frac{∣ \nabla _{rad, ϕ} ∣}{∣ \nabla _{ad, ϕ} ∣} - 1)_{background} = \frac{1/ ( 4 - β )}{( γ - 1 ) / ( γ )} - 1,

(\frac{∣ \nabla _{rad, ϕ} ∣}{∣ \nabla _{ad, ϕ} ∣} - 1)_{background} = \frac{1/ ( 4 - β )}{( γ - 1 ) / ( γ )} - 1,

κ (\overset{ˉ}{T} (r, ϕ)) = κ_{0} \overset{ˉ}{T}^{β} .

κ (\overset{ˉ}{T} (r, ϕ)) = κ_{0} \overset{ˉ}{T}^{β} .

κ (r) = 1.11 (\frac{r}{a _{p}})^{- δ} c m^{2} g^{- 1},

κ (r) = 1.11 (\frac{r}{a _{p}})^{- δ} c m^{2} g^{- 1},

T^{n + 1} = η_{1} + η_{2} E_{R}^{n + 1},

T^{n + 1} = η_{1} + η_{2} E_{R}^{n + 1},

η_{1} = \frac{T ^{n} + 12Δ t \frac{κ _{P}}{c _{V}} σ ( T ^{n} ) ^{4} + Δ t \frac{Q _{ϵ - indep}}{ρ c _{V}}}{1 + 16Δ t \frac{κ _{P}}{c _{V}} σ ( T ^{n} ) ^{3} + Δ t ( γ - 1 ) \nabla \cdot \@vec v},

η_{1} = \frac{T ^{n} + 12Δ t \frac{κ _{P}}{c _{V}} σ ( T ^{n} ) ^{4} + Δ t \frac{Q _{ϵ - indep}}{ρ c _{V}}}{1 + 16Δ t \frac{κ _{P}}{c _{V}} σ ( T ^{n} ) ^{3} + Δ t ( γ - 1 ) \nabla \cdot \@vec v},

η_{2} = \frac{Δ t \frac{κ _{P}}{c _{V}} c}{1 + 16Δ t \frac{κ _{P}}{c _{V}} σ ( T ^{n} ) ^{3} + Δ t ( γ - 1 ) \nabla \cdot \@vec v} .

η_{2} = \frac{Δ t \frac{κ _{P}}{c _{V}} c}{1 + 16Δ t \frac{κ _{P}}{c _{V}} σ ( T ^{n} ) ^{3} + Δ t ( γ - 1 ) \nabla \cdot \@vec v} .

(\@vec a \cdot \nabla) \@vec a = \frac{1}{2} \nabla (\@vec a \cdot \@vec a) + (\nabla \times \@vec a) \times \@vec a,

(\@vec a \cdot \nabla) \@vec a = \frac{1}{2} \nabla (\@vec a \cdot \@vec a) + (\nabla \times \@vec a) \times \@vec a,

\nabla \times (\@vec a \times \@vec b) = \@vec a (\nabla \cdot \@vec b) - \@vec b (\nabla \cdot \@vec a) + (\@vec b \cdot \nabla) \@vec a - (\@vec a \cdot \nabla) \@vec b,

\nabla \times (\@vec a \times \@vec b) = \@vec a (\nabla \cdot \@vec b) - \@vec b (\nabla \cdot \@vec a) + (\@vec b \cdot \nabla) \@vec a - (\@vec a \cdot \nabla) \@vec b,

\nabla \cdot (\nabla \times \@vec a) = 0,

\nabla \cdot (\nabla \times \@vec a) = 0,

\nabla \times (\nabla f) = \@vec 0,

\nabla \times (\nabla f) = \@vec 0,

\nabla \times [(\@vec v \cdot \nabla) \@vec v] = \nabla \times [\frac{1}{2} \nabla (\@vec v \cdot \@vec v) + (\nabla \times \@vec v) \times \@vec v] = \nabla \times [\@vec ω_{r} \times \@vec v],

\nabla \times [(\@vec v \cdot \nabla) \@vec v] = \nabla \times [\frac{1}{2} \nabla (\@vec v \cdot \@vec v) + (\nabla \times \@vec v) \times \@vec v] = \nabla \times [\@vec ω_{r} \times \@vec v],

\nabla \times [(\@vec v \cdot \nabla) \@vec v] = \@vec ω_{r} (\nabla \cdot \@vec v) + (\@vec v \cdot \nabla) \@vec ω_{r} - (\@vec ω_{r} \cdot \nabla) \@vec v,

\nabla \times [(\@vec v \cdot \nabla) \@vec v] = \@vec ω_{r} (\nabla \cdot \@vec v) + (\@vec v \cdot \nabla) \@vec ω_{r} - (\@vec ω_{r} \cdot \nabla) \@vec v,

\nabla \times (\frac{\nabla P}{ρ}) = - \frac{1}{ρ ^{2}} \nabla ρ \times \nabla P,

\nabla \times (\frac{\nabla P}{ρ}) = - \frac{1}{ρ ^{2}} \nabla ρ \times \nabla P,

\nabla \times (2 \@vec Ω \times \@vec v) = 2 [\@vec Ω (\nabla \cdot \@vec v) - (\@vec Ω \cdot \nabla) \@vec v],

\nabla \times (2 \@vec Ω \times \@vec v) = 2 [\@vec Ω (\nabla \cdot \@vec v) - (\@vec Ω \cdot \nabla) \@vec v],

\frac{\partial \@vec ω _{r}}{\partial t} + (\@vec v \cdot \nabla) \@vec ω_{r} = \frac{D \@vec ω _{r}}{D t} = (\@vec ω_{a} \cdot \nabla) \@vec v - \@vec ω_{a} (\nabla \cdot \@vec v) + \frac{\nabla ρ \times \nabla P}{ρ ^{2}},

\frac{\partial \@vec ω _{r}}{\partial t} + (\@vec v \cdot \nabla) \@vec ω_{r} = \frac{D \@vec ω _{r}}{D t} = (\@vec ω_{a} \cdot \nabla) \@vec v - \@vec ω_{a} (\nabla \cdot \@vec v) + \frac{\nabla ρ \times \nabla P}{ρ ^{2}},

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

11institutetext: Institute of Astronomy, Charles University in Prague, V Holešovičkách 2, CZ–18000 Prague 8, Czech Republic

11email: [email protected] 22institutetext: Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, 22100 Lund, Sweden

Oscillatory migration of accreting protoplanets

driven by a 3D distortion of the gas flow

O. Chrenko 11

M. Lambrechts 22

(Received ??; accepted ??)

Abstract

*Context. * The dynamics of a low-mass protoplanet accreting solids is influenced by the heating torque, which was found to suppress inward migration in protoplanetary disks with constant opacities.

*Aims. * We investigate the differences of the heating torque between disks with constant and temperature-dependent opacities.

*Methods. * Interactions of a super-Earth-sized protoplanet with the gas disk are explored using 3D radiation hydrodynamic simulations.

*Results. * Accretion heating of the protoplanet creates a hot underdense region in the surrounding gas, leading to misalignment of the local density and pressure gradients. As a result, the 3D gas flow is perturbed and some of the streamlines form a retrograde spiral rising above the protoplanet. In the constant-opacity disk, the perturbed flow reaches a steady state and the underdense gas responsible for the heating torque remains distributed in accordance with previous studies. If the opacity is non-uniform, however, the differences in the disk structure can lead to more vigorous streamline distortion and eventually to a flow instability. The underdense gas develops a one-sided asymmetry which circulates around the protoplanet in a retrograde fashion. The heating torque thus strongly oscillates in time and does not on average counteract inward migration.

*Conclusions. * The torque variations make the radial drift of the protoplanet oscillatory, consisting of short intervals of alternating rapid inward and outward migration. We speculate that transitions between the positive and oscillatory heating torque may occur in specific disk regions susceptible to vertical convection, resulting in the convergent migration of multiple planetary embryos.

Key Words.:

Hydrodynamics – Planets and satellites: formation – Planet-disk interactions – Protoplanetary disks

1 Introduction

Migration of protoplanets embedded in their natal gas disks is a key evolutionary step in formation of each planetary system. Low-mass protoplanets, incapable of gap opening (Crida et al., 2006), undergo Type I migration under the influence of the gravitational torques exerted by the Lindblad spiral wakes (Goldreich & Tremaine, 1979; Ward, 1986) and by the gas in their corotation region (Ward, 1991; Masset, 2002; Tanaka et al., 2002; Paardekooper & Mellema, 2006; Baruteau & Masset, 2008; Paardekooper & Mellema, 2008; Masset & Casoli, 2009; Paardekooper & Papaloizou, 2009; Baruteau et al., 2011). Type I migration depends in a complicated way on the disk structure and thermophysics (e.g. Kley & Crida, 2008; Paardekooper et al., 2010, 2011; Lega et al., 2015). Its detailed understanding is therefore essential for the creation of realistic population synthesis models (e.g. Coleman & Nelson, 2016).

It has been recently discovered that low-mass protoplanets evolving in radiative disks are subject to thermal torques (Lega et al., 2014; Benítez-Llambay et al., 2015; Masset & Velasco Romero, 2017; Masset, 2017) related to the thermal perturbations induced by a protoplanet in its vicinity. If the protoplanet itself is cold and non-luminous (Lega et al., 2014), the gas which arrives into its potential well becomes heated mostly as a result of compression (by means of the thermodynamic ‘ $P\mathrm{d}V$ ’ term). The arising temperature excess becomes smoothed out by the radiative transfer, so when the gas leaves the high pressure region it lacks some of its internal energy – compared to the state before the compression – and therefore becomes colder and overdense. Two overdense lobes appear along the streamlines outflowing from the Hill sphere and their asymmetry makes the total torque felt by the protoplanet more negative, enhancing the inward migration. The process is known as the cold-finger effect (Lega et al., 2014).

In the opposite limit, the protoplanet is hot, as a result of the solid material deposition during its formation (e.g. by pebble accretion; Ormel & Klahr, 2010; Lambrechts & Johansen, 2012). In such a case, the luminous protoplanet acts as a local heat source for the surrounding gas. Once the gas is heated, it becomes underdense compared to the situation without accretion heating. Benítez-Llambay et al. (2015) performed 3D radiation hydrodynamic simulations with the assumption of the constant disk opacity and found that the hot protoplanet on a fixed circular orbit creates two underdense lobes of gas, again associated with the outflow from the Hill sphere. The rear lobe (positioned behind the protoplanet outwards from its orbit) is dominant, therefore there is an overabundance of gas ahead of the protoplanet and the resulting torque becomes more positive, supporting outward migration. The positive enhancement was named the heating torque. It was proposed to be an additional mechanism (along with other posibilites, see e.g. Rafikov, 2002; Paardekooper & Mellema, 2006; Morbidelli et al., 2008; Li et al., 2009; Yu et al., 2010; Kretke & Lin, 2012; Bitsch et al., 2013; Fung & Chiang, 2017; Brasser et al., 2018; Miranda & Lai, 2018; McNally et al., 2019) capable of preventing the destruction of terrestrial-sized planetary embryos by an overly efficient inward migration (Korycansky & Pollack, 1993; Ward, 1997; Tanaka et al., 2002).

Moreover, the heating torque has important dynamical consequences for migrating protoplanets (Brož et al., 2018; Chrenko et al., 2018) because it can excite orbital eccentricities and inclinations by means of the hot-trail effect (Eklund & Masset, 2017; Chrenko et al., 2017), which counteracts the otherwise efficient eccentricity and inclination damping by waves (Tanaka & Ward, 2004; Cresswell et al., 2007).

Nevertheless, the heating torque has not been extensively studied in 3D radiative disks with non-uniform opacities. However, realistic opacity functions (Bell & Lin, 1994; Semenov et al., 2003; Zhu et al., 2012) are of a great significance for the disk structure and planet migration (e.g. Kretke & Lin, 2012; Bitsch et al., 2013) and this work will show that the heating torque is affected as well.

In this paper, we reinvestigate the thermal torques acting on a low-mass protoplanet, with a special emphasis on the heating torque in a disk with non-uniform opacity. We examine the streamlines near the protoplanet and point out the importance of their 3D distortion for redistribution of the hot underdense gas responsible for the heating torque.

2 Model

We consider a protoplanetary system consisting of a central protostar surrounded by a disk of coupled gas and dust in which a single protoplanet is embedded. The fluid part of the disk model accounts only for the gas, assuming the dust is a passive tracer that acts as the main contributor to the material opacity.

2.1 Governing equations

The disk is described using Eulerian hydrodynamics on a staggered spherical mesh centred on the protostar. The spherical coordinates consist of the radial distance $r$ , azimuthal angle $\theta$ and colatitude $\phi$ . Our model is built on top of the hydrodynamic module of fargo3d111Public version of the code is available at http://fargo.in2p3.fr/. (Benítez-Llambay & Masset, 2016), which solves the equations of continuity and momentum of a fluid

[TABLE]

where the individual quantities stand for the volume density $\rho$ , time $t$ , flow velocity vector $\@vec{v}$ (with the radial, azimuthal and vertical components $v_{r}$ , $v_{\theta}$ and $v_{\phi}$ ), pressure $P$ , gravitational potential $\Phi$ of the protostar and the protoplanet, viscous stress tensor $\@tens{T}$ (Benítez-Llambay & Masset, 2016), radius vector $\@vec{r}$ and angular velocity vector $\@vec{\Omega}$ which is non-zero because we work in the reference frame corotating with the protoplanet.

To account for the energy production, dissipation and transport in the disk, we implement the two-temperature energy equations for the gas and the radiation field according to the formulation of Bitsch et al. (2013):

[TABLE]

where $E_{\mathrm{R}}$ is the radiative energy, $\@vec{F}$ the radiation flux, $\kappa_{\mathrm{P}}$ the Planck opacity, $\sigma$ the Stefan-Boltzmann constant, $T$ the gas temperature, $c$ the speed of light, $\epsilon$ the internal energy of the gas, $Q_{\mathrm{visc}}$ the viscous heating term (Mihalas & Weibel Mihalas, 1984), $Q_{\mathrm{art}}$ describes the heating due to the artificial viscosity (Stone & Norman, 1992) and $Q_{\mathrm{acc}}$ the heat released when the protoplanet is accreting. Stellar irradiation is neglected in this paper for simplicity although our code is capable of including it as well.

The state equation of the ideal gas is used:

[TABLE]

where $\gamma$ is the adiabatic index and $c_{V}$ is the specific heat at constant volume, which can be expressed as $c_{V}=R/(\mu(\gamma-1))$ , where $R$ is the universal gas constant and $\mu$ is the mean molecular weight.

The flux-limited diffusion approximation (FLD; Levermore & Pomraning, 1981) is adopted to obtain a closure relation for the radiation flux

[TABLE]

where $\kappa_{\mathrm{R}}$ is the Rosseland opacity and $\lambda_{\mathrm{lim}}$ is the flux limiter of Kley (1989). For the opacities, we assume that the Planck and Rosseland means are similar enough to be replaced with a single value $\kappa$ . This is a valid assumption in the cold regions of protoplanetary disks that we aim to study (Bitsch et al., 2013). The detailed opacity law will be specified later in Sect. 2.3.

The accretion luminosity of the protoplanet is given by

[TABLE]

where $G$ is the gravitational constant, $M_{\mathrm{p}}$ is the mass of the protoplanet, $\dot{M}_{\mathrm{p}}$ is its accretion rate, and $R_{\mathrm{p}}$ is the protoplanet’s radius. In writing the second equality, we introduce the mass doubling time $\tau=M_{\mathrm{p}}/\dot{M}_{\mathrm{p}}$ , which is a free parameter that controls the accretion rate in our model (see Sect. 2.2).

We assume that the radiation flux from the protoplanet is completely absorbed by the optically thick gas in the eight grid cells enclosing the protoplanet (Benítez-Llambay et al., 2015; Eklund & Masset, 2017; Lambrechts & Lega, 2017). The accretion heating, which is non-zero only within these cells, is then simply

[TABLE]

where $V$ is the total volume of the heated cells.

The disk evolves in the combined gravitational potential of the central protostar and the embedded protoplanet:

[TABLE]

where $M_{\star}$ is the mass of the protostar and $d$ is the cell-protoplanet distance. The planetary potential is smoothed to avoid numerical divergence at the protoplanet location ( $d=0$ ) using the tapering cubic-spline function $f_{\mathrm{sm}}$ of Klahr & Kley (2006):

[TABLE]

where the smoothing length $r_{\mathrm{sm}}$ is a fraction of the protoplanet’s Hill sphere radius (see Sect. 2.2). Orbital evolution of the protoplanet is tracked using the ias15 integrator (Rein & Spiegel, 2015) from the rebound222Public version of the code is available at https://rebound.readthedocs.io/en/latest/. package (Rein & Liu, 2012), which we interfaced with fargo3d.

Although fargo3d is an explicit hydrodynamic solver, implementing Eqs. (3) and (4) in an explicit form would lead to a very restrictive Courant condition on the largest allowed time step length. To avoid such a time step limitation, it is advantageous to solve the energy equations in an implicit form. We thus follow the discretisation and linearisation proposed by Bitsch et al. (2013), with a minor modification introduced in Appendix A.

The solution of the implicit problem is obtained iteratively, by minimizing the residual $\@vec{r}\equiv\mathbf{A}\@vec{x}-\@vec{b}$ , where $\mathbf{A}$ is the matrix of the linear system, $\@vec{x}$ is the solution vector and $\@vec{b}$ is the right-hand side vector. Our iterative solver uses the improved bi-conjugate stabilised method (IBiCGStab; Yang & Brent, 2002) with the Jacobi preconditioning. The convergence criterion is $\|{\@vec{r}}\|/\|{\@vec{b}}\|<10^{-4}$ , where the norms are calculated in the $\mathrm{L}_{2}$ space.

2.2 Initial conditions and parameters

Initial conditions describing the disk are adopted from Kley et al. (2009) and Lega et al. (2014). The surface density follows the power-law function $\Sigma=484(r/(1\,\mathrm{au}))^{-0.5}\,\mathrm{g\,cm^{-2}}$ . The initial $\Sigma$ is converted into the initial $\rho$ assuming the disk is vertically isothermal. The velocity field is given by balancing the gravitational acceleration from the protostar, the acceleration due to the pressure gradient and the centrifugal acceleration. We assign constant kinematic viscosity $\nu=10^{15}\,\mathrm{cm^{2}\,s^{-1}}$ to the gas to mimic the angular momentum transport in the disk driven by physical effects outside the scope of this study (e.g. the turbulent eddy viscosity). The initial aspect ratio is $h=H/r=0.05$ , where $H$ is the pressure scale height. The disk is therefore initially non-flaring but we point our that $h$ evolves during the simulations. The adiabatic index is $\gamma=1.43$ and the mean molecular weight is $\mu=2.3\,\mathrm{g\,mol^{-1}}$ , corresponding to the solar mixture of the hydrogen and helium.

We will assume a single embedded super-Earth with the mass $M_{\mathrm{p}}=3\,M_{\oplus}$ , orbiting on a circular orbit at the Jupiter’s heliocentric distance $a_{\mathrm{p}}\equiv a_{\mathrm{J}}=5.2\,\mathrm{au}$ . Unless otherwise specified, the orbit is held fixed in our simulations, i.e. the protoplanet is not allowed to migrate and the torque we measure is the static torque. The simulations are performed in the reference frame corotating with the protoplanet at the Keplerian angular frequency $\Omega_{\mathrm{p}}\equiv\Omega_{\mathrm{J}}=\sqrt{GM_{\star}/a_{\mathrm{J}}^{3}}$ and the smoothing length for the tapering function of the gravitational potential $f_{\mathrm{sm}}$ is $r_{\mathrm{sm}}=0.5R_{\mathrm{H}}$ , where $R_{\mathrm{H}}=a_{\mathrm{p}}(M_{\mathrm{p}}/(3M_{\star}))^{1/3}$ is the Hill sphere radius. We will study both cases without and with planetary accretion, corresponding to the cold- and hot-protoplanet limit, respectively. In the simulations with the accretion, we assume the mass doubling time $\tau=100\,\mathrm{kyr}$ , which is a value within the range of the expected pebble accretion rates (Lambrechts & Johansen, 2014; Chrenko et al., 2017). The resulting luminosity of the protoplanet is $L\simeq 4.2\times 10^{27}\,\mathrm{erg}\,\mathrm{s}^{-1}$ .

2.3 Opacities

Using the initial disk parameters, we setup two fiducial disk models that only differ in the material opacity function. The opacity is either constant or follows a slightly modified (see below) prescription of Bell & Lin (1994). To distinguish between the models, we use the abbreviations $\kappa_{\mathrm{const}}$ -disk and $\kappa_{\mathrm{BL}}$ -disk, respectively. The exact value of the opacity in the $\kappa_{\mathrm{const}}$ -disk is tuned to be the same as in the $\kappa_{\mathrm{BL}}$ -disk at the location of the protoplanet $a_{\mathrm{p}}$ , leading to $\kappa_{\mathrm{const}}\equiv\kappa_{\mathrm{BL}}(r=a_{\mathrm{p}})=1.11\,\mathrm{cm^{2}\,g^{-1}}$ .

The unmodified opacity from Bell & Lin (1994), which we denote $\kappa_{\mathrm{BL}}^{\mathrm{full}}(\rho(r,\theta,\phi),T(r,\theta,\phi))$ , is a fitting law that sets $\kappa$ as a function of the local gas density and temperature. The table spans several regimes corresponding to the presence (or absence) of dust or molecular species dominant in protoplanetary disks. Using $\kappa_{\mathrm{BL}}^{\mathrm{full}}$ in our simulations with the accretion heating of the protoplanet could cause strong local opacity gradients, because we expect the temperature perturbations to reach $\sim$$10^{1}\,\mathrm{K}$ at distances of one cell size from the centre of the protoplanet.

In practice, we apply the Bell & Lin’s opacity law in a simplified way, using $\kappa_{\mathrm{BL}}(\bar{\rho}(r,\phi),\bar{T}(r,\phi))$ , where the bared quantities are azimuthally computed arithmetic means and the dependence on the $\theta$ -coordinate is therefore dropped. This helps us to distinguish the effects caused by the global structure of the disk from those related to the local $\kappa$ - $T$ - $\rho$ feedback. To justify our approach, we verify in Appendix B that the unmodified opacity table $\kappa_{\mathrm{BL}}^{\mathrm{full}}$ of Bell & Lin (1994) does not change our conclusions.

The fiducial $\kappa_{\mathrm{const}}$ - and $\kappa_{\mathrm{BL}}$ -disks are discussed throughout the majority of the paper, with the exception of Sect. 3.6 where we study the dependence of the heating torque on the opacity gradient within the disk.

2.4 Grid resolution and boundary conditions

Migration of low-mass protoplanets critically depends on the grid resolution, we therefore combine the well-established disk extent and resolution from Lega et al. (2014) (in the azimuthal and vertical direction) and Eklund & Masset (2017) (in the radial direction). The disk radially stretches from $r_{\mathrm{min}}=3.12\,\mathrm{au}$ to $r_{\mathrm{max}}=7.28\,\mathrm{au}$ and is resolved by 512 rings. We prevent any vertical motions of the protoplanet in our simulations by assuming that the solution is symmetric with respect to the midplane. One of the disk boundaries in the colatitude is therefore located at the midplane ( $\phi=\pi/2$ ); the vertical extent above the midplane is $7^{\circ}$ . The colatitude is resolved by 64 zones. In the azimuth, we use only 1 sector in our preparatory simulations without the protoplanet and 1382 sectors in our simulations with the embedded protoplanet. The resulting local resolution is 8 cells per $R_{\mathrm{H}}$ in the $r$ - and $\phi$ -directions and 3 cells per $R_{\mathrm{H}}$ in the $\theta$ -direction (see also Appendix B.2 for a simple resolution test).

The azimuthal boundary conditions are periodic for all primitive quantities. The radial boundary conditions are symmetric for $\rho$ , $\epsilon$ , and $v_{\phi}$ and reflecting for $v_{r}$ . $E_{\mathrm{R}}$ is set to a zero gradient and $v_{\theta}$ is extrapolated using the same radial dependence as for the Keplerian rotation velocity. The boundary conditions in colatitude are symmetric for $\rho$ , $\epsilon$ , $v_{r}$ and $v_{\theta}$ , and reflecting for $v_{\phi}$ . $E_{\mathrm{R}}$ is symmetrised at the midplane boundary and set to $a_{\mathrm{R}}T_{0}^{4}$ at the remaining boundary in the colatitude, where $a_{\mathrm{R}}$ is the radiation constant and $T_{0}\equiv 5\,\mathrm{K}$ is the ambient temperature that allows for vertical radiative cooling of the disk. Additionally, wave-damping conditions of de Val-Borro et al. (2006) are imposed near the radial edges and also near the disk surface.

3 Simulations

3.1 Equilibrium disks

Since we use a non-isothermal equation of state and we also account for the energy production and transfer, the parametric setup of the disk, which was discussed so far, is not stationary. Therefore, before the simulations with the embedded protoplanet are conducted, we let both disks ( $\kappa_{\mathrm{const}}$ -disk and $\kappa_{\mathrm{BL}}$ -disk) numerically evolve over the time scale $100\,P_{\mathrm{orb}}$ , where $P_{\mathrm{orb}}=2\pi/\Omega_{\mathrm{J}}$ .

The equilibrium state after the relaxation is characterized by Fig. 1. We plot the radial profiles of the midplane temperature $T$ , opacity $\kappa$ , entropy $S=P/\rho^{\gamma}$ and aspect ratio $h=H/r$ , where the pressure scale height is determined as $H=c_{\mathrm{s}}/\Omega_{\mathrm{K}}$ and the sound speed as $c_{\mathrm{s}}=\sqrt{\gamma P/\rho}$ .

The constant opacity of the $\kappa_{\mathrm{const}}$ -disk makes all the remaining radial dependencies rather shallow. The aspect ratio is almost flat, only slightly increasing with the radial distance. In the $\kappa_{\mathrm{BL}}$ -disk, on the other hand, the opacity has a peak near $3.5\,\mathrm{au}$ , where the water ice evaporates for the given setup, and decreases at larger radii. Therefore, the efficiency of the disk cooling increases and simultaneously, the efficiency of the viscous heating diminishes with the dropping relative velocity of the shearing layers. Consequently, the aspect ratio radially decreases as the energy budget is not sufficient to puff up the disk. In such a disk, $T$ and $S$ profiles radially decrease more steeply compared to the $\kappa_{\mathrm{const}}$ -disk.

3.2 Torque evolution

Starting from the relaxed state of the disks, we copy the hydrodynamic quantities in the azimuthal direction to expand the resolution from a single sector to the desired 1382 sectors. The protoplanet is inserted and we simulate $30\,P_{\mathrm{orb}}$ of evolution while neglecting any accretion and accretion heating of the protoplanet. The aim of this part of the simulation is to let the disk adjust to the presence of a gravitational perturber and to acquire a converged value of the disk torque in the absence of the accretion heating, corresponding to the cold-protoplanet limit. For the $\kappa_{\mathrm{BL}}$ -disk, this part of the simulation is similar to the experiments in Lega et al. (2014). Subsequently, we continue the simulation for another $30\,P_{\mathrm{orb}}$ during which we let the planetary mass grow while releasing the accretion heat into the gas disk according to Eqs. (7) and (8).

Fig. 2 shows the temporal evolution of the torque exerted on the protoplanet by the $\kappa_{\mathrm{const}}$ -disk and $\kappa_{\mathrm{BL}}$ -disk. During the phase without accretion heating, the torque converges to a stationary value during $10\,P_{\mathrm{orb}}$ . The torque value in the $\kappa_{\mathrm{BL}}$ -disk is more positive compared to the $\kappa_{\mathrm{const}}$ -disk, which is because the steeper radial decline of the entropy in the $\kappa_{\mathrm{BL}}$ -disk enhances the positive entropy-driven part of the corotation torque (Paardekooper & Mellema, 2006; Baruteau & Masset, 2008).

When the accretion heating is activated in the $\kappa_{\mathrm{const}}$ -disk, a positive contribution is added to the torque, in accordance with the results of Benítez-Llambay et al. (2015). The torque slightly oscillates at first, but the amplitude of the oscillations decreases in time and becomes negligible at $t=50\,P_{\mathrm{orb}}$ .

On contrary, when the accretion heating is activated in the $\kappa_{\mathrm{BL}}$ -disk, the outcome of the heating torque becomes very different. Strong oscillations of the disk torque are excited almost immediately and they do not vanish in time; instead, their amplitude remains the same. The arithmetic mean of the torque measured in the time interval $30$ – $60\,P_{\mathrm{orb}}$ is $\bar{\Gamma}\simeq-6.3\times 10^{-6}\,a_{\mathrm{J}}^{2}\,\Omega_{\mathrm{J}}^{2}$ , implying that the torque is more negative compared to the situation without accretion heating (the mean torque in the interval $20$ – $30\,P_{\mathrm{orb}}$ is $\bar{\Gamma}\simeq-1.6\times 10^{-6}\,a_{\mathrm{J}}^{2}\,\Omega_{\mathrm{J}}^{2}$ ). The amplitude of the variations with respect to the mean value is $\simeq$$\pm 3.9\times 10^{-5}\,a_{\mathrm{J}}^{2}\,\Omega_{\mathrm{J}}^{2}$ and the oscillation period is $\simeq$$2.1\,P_{\mathrm{orb}}$ .

The torque oscillations are unexpected and we thus focus throughout the rest of the paper on finding the physical mechanism that excites them. The occurrence of oscillations suggests that the gas distribution around the protoplanet is changing during the simulation, as we can demonstrate using the integral expression for the disk torque

[TABLE]

where $\@vec{r}_{\mathrm{p}}$ is the protoplanet’s radius vector, $\@vec{F}_{\mathrm{g}}$ is the gravitational force of a disk element, the vertical component of the cross product is considered and we integrate over the disk volume. Only a non-zero azimuthal component of $\@vec{F}_{\mathrm{g}}$ can lead to a non-vanishing cross product in the integral, thus any oscillations of $\Gamma$ must be related to a variation of $F_{\mathrm{g},\theta}$ . In other words, there must be an azimuthal asymmetry in the gas distribution with respect to the protoplanet for the torque to be non-zero and only a temporal redistribution of the asymmetry can cause a torque oscillation.

Our strategy throughout the rest of Sect. 3 is the following: First, we focus on the $\kappa_{\mathrm{const}}$ -disk simulation in Sect. 3.3. Although similar simulations were analysed by Benítez-Llambay et al. (2015), our aim is to focus on the 3D gas flow that has not been described yet. Our findings are then expanded for the $\kappa_{\mathrm{BL}}$ -disk simulation in Sect. 3.4 where we relate the gas redistribution to the oscillatory behaviour of the torque. Sect. 3.5 is devoted to identifying the key physical mechanisms that affect the gas flow. In Sect. 3.6, we vary the disk opacity gradient and we study its impact on the torque oscillation. Finally, we relax the assumption of a fixed orbit and explore how the protoplanet migrates in Sect. 3.7.

3.3 Steady state of the heated gas

In this section, the $\kappa_{\mathrm{const}}$ -disk simulation is analysed.

3.3.1 Midplane

Fig. 3 compares the midplane density and temperature distribution around the cold and hot protoplanet. We display the state of the simulation at $t=30$ (i.e. at the final stage of the phase without accretion heating) and $60\,P_{\mathrm{orb}}$ (i.e. at the final stage of the phase with accretion heating). The latter represents the steady state of gas around the accreting protoplanet and we verified that such a distribution is achieved early (at $t\simeq 31\,P_{\mathrm{orb}}$ ) and does not greatly evolve since then.

For the non-luminous protoplanet, the gas state is in agreement with Lega et al. (2014) (see their Fig. 10 for a comparison). The density distribution is not spherically symmetric with respect to the protoplanet but there are two patches of slightly overdense gas along the outflow from the Hill sphere known as the cold fingers (explained in Sect. 1). The temperature drop inside the fingers is clearly apparent from the top right panel of Fig. 3.

When the protoplanet becomes luminous, the gas distribution is modified and the heating torque arises. Following Benítez-Llambay et al. (2015) and Masset (2017), one can imagine the response of the gas to the heating from the protoplanet as follows: First, an underdense disturbance appears close to the protoplanet, with a characteristic scale length given by the linear perturbation model of Masset (2017):

[TABLE]

where $\chi$ is the thermal diffusivity and $q$ is a dimensionless measure of the disk shear ( $q=3/2$ for a Keplerian disk). Second, the low-density gas is distorted by the shear motions. The rotation of the inner disk with respect to the protoplanet is faster and the low-density gas thus propagates ahead of the protoplanet. The motion of the outer disk lags behind the protoplanet and so does the hot perturbation. As a result, two hot lobes with decreased density are formed along the streamlines outflowing from the Hill sphere, as described by Benítez-Llambay et al. (2015). The size of the lobes is inherently asymmetric because the corotation between the protoplanet and the sub-Keplerian gas is radially shifted inwards, therefore the outer rear lobe is usually larger and the heating torque should be positive.

Such an advection-diffusion interplay is indeed observed in Fig. 3, where we find the typical two-lobed distribution of hot underdense gas around the luminous protoplanet and the positive boost of the total torque (Fig. 2) confirms that the outer lobe is slightly more pronounced. The bottom right panel reveals the magnitude and spatial extent of the temperature excess, as well as its skewed shape in the direction of the disk shear.

However, we make a new observation here concerning the streamlines of the flow that are overlaid in the temperature maps. It is obvious that the hot perturbation significantly changes the topology of the flow with respect to the cold-protoplanet case. U-turn streamlines no longer appear in the depicted part of the disk, the direction of the circulating streamlines changes as they pass the protoplanet and a new set of spiral-like retrograde streamlines appears.

3.3.2 Vertical plane

It is known that vertical motions play an important role for the structure of circumplanetary envelopes (e.g. Tanigawa et al., 2012; Fung et al., 2015; Ormel et al., 2015; Cimerman et al., 2017; Lambrechts & Lega, 2017; Kurokawa & Tanigawa, 2018; Popovas et al., 2019). Since previous studies of the heating torque did not investigate the vertical perturbations, we thus do so here. Fig. 4 shows the gas temperature and the velocity field in the vertical plane intersecting the protoplanet’s location. When the protoplanet is cold, there is a vertical stream of gas descending towards the protoplanet and escaping as an outflow through the midplane, in accordance with e.g. Lambrechts & Lega (2017).

In the presence of accretion heating, however, the direction of the gas flow above the protoplanet is reverted. It forms an outflowing column while the midplane flow becomes directed towards the protoplanet. Two overturning cells appear on each side of the vertical column (although it is important to point out that no such cells are apparent in the full 3D flow which is discussed later). We notice that the hot perturbation is not spherically symmetric but rather elongated in the direction of the column, indicating that the envelope is not hydrostatic.

3.3.3 2D and 3D streamline topology

So far, we described two new findings that were not incorporated in the existing descriptions of the heating torque (Benítez-Llambay et al., 2015; Masset, 2017): The distortion of the streamline topology and the reversal of the vertical motions. The flow direction is directly linked to the heating torque because it determines the redistribution of the hot gas by advection and thus contributes to the shape of the underdense regions. Therefore we focus on the streamline topology in this section.

The streamlines are calculated using the explicit first order Euler integrator and the trilinear interpolation of the velocity field. The interpolation allows us to obtain the velocity vector at an arbitrary location within the spherical grid. The size of the integration step is chosen so that the length integrated during a single propagation does not exceed $0.1$ of the shortest cell dimension.

We construct 2D and 3D projections of the streamline topologies. For the 2D projections, the streamlines are generated exactly at the midplane where $v_{\phi}=0$ and although they provide a useful visualisation, we emphasise that they – by construction – carry zero information about the adjacent vertical flows. For the 3D projections, the streamlines are generated slightly above the midplane (at $\phi=\pi/2-0.005\,\mathrm{rad}$ ) to take into account non-zero $v_{\phi}$ .

Fig. 5 shows 2D and 3D streamlines in the $\kappa_{\mathrm{const}}$ -disk, again for the exact same simulation stages as in Figs. 3 and 4. In the plots, we distinguish the following types of streamline: First, there are circulating streamlines that do not cross the protoplanet’s corotation. To imagine the direction of the relative motion, we point out that the gas on inner circulating streamlines moves faster than the protoplanet whereas the gas on outer circulating streamlines lags behind the protoplanet. Second, there are horseshoe streamlines that make a single U-turn and cross the corotation once at each side of the protoplanet. Such streamlines ahead of the protoplanet’s orbital motion form the front horseshoe region, those located behind the protoplanet belong to the rear horseshoe region. To outline the separatrices between the horseshoe and circulating regions, we highlight the critical inner and outer circulating streamlines that are located closest to the protoplanet. Finally, some of our plots contain streamlines that do not fall in any of the aforementioned categories.

When the protoplanet is non-luminous, the 2D midplane streamlines in Fig. 5 do not exhibit any unexpected features. The stagnation point (X-point) of the flow is located within the Hill sphere, which is also intersected by both horseshoe and circulating streamlines. In 3D, we notice that upon making their U-turn, the horseshoe streamlines vertically descend towards the midplane, as already pointed out by Fung et al. (2015) or Lambrechts & Lega (2017). The descend is also exhibited by some of the circulating streamlines closest to the protoplanet.

For the hot protoplanet, we now obtain a clear picture of the streamline distortion. In the midplane, the following changes appear:

•

Circulating streamlines cross a smaller portion of the Hill sphere. When passing the protoplanet, they are bent towards it (unlike near the non-luminous protoplanet where they are rather deflected away), which is especially apparent for the critical circulating streamlines.

•

The classical horseshoe streamlines detach from the Hill sphere and make their U-turns at greater azimuthal separations.

•

Part of the streamlines originating downstream the horseshoe regions is captured inside the Hill sphere where it rotates around the protoplanet in a retrograde fashion.

In 3D, the distortion has the following additional features:

•

The ‘captured’ streamlines are uplifted and form a spiral-like vertical column, outflowing and escaping from the Hill sphere.

•

When the circulating streamlines pass the protoplanet and become perturbed, they are also uplifted. This behaviour is exactly the opposite to the cold-protoplanet situation.

3.4 Instability of the heated gas

We now return to the $\kappa_{\mathrm{BL}}$ -disk simulation for which we discovered the strong oscillations of the heating torque (Fig. 2).

3.4.1 Evolving underdense lobes

Investigating the evolution of the gas density, we find that the position and size of the underdense lobes never becomes stationary, as shown in Fig. 6 (see also the online movie 1). The figure shows a sequence of snapshots corresponding to $t=31.3$ , $31.75$ , $32.20$ and $32.65\,P_{\mathrm{orb}}$ .

The first panel depicts the state when the total torque reaches its first minimum during the beginning of the oscillatory phase. There is a dominant underdense lobe located ahead of the protoplanet while the rear lobe almost disappears. Such a distribution can be easily related to the strong negative torque: The excess of the gas mass behind the protoplanet (and the paucity of mass ahead of it) leads to an azimuthal pull acting against the orbital motion, imposing a negative torque.

When the torque is reversing from negative to positive (second panel), both lobes are similarly pronounced. The rear one seems to be located closer to the protoplanet. The third panel corresponds to the torque maximum. The rear lobe is dominant and thus the overabundance of the gas ahead of the protoplanet makes the torque positive. The final panel shows the state when the oscillating torque is approximately half-way from positive to negative. The gas distribution indeed looks like a counterpart to the second panel since both lobes are again apparent but the front one is now closer to the protoplanet.

The gas redistribution is clearly related to the topology of the streamlines and to the position of the spiral-like flow. We notice that the centre of the captured streamlines undergoes retrograde (‘clockwise’) rotation around the protoplanet. One underdense lobe is always associated with this rotating flow. Apparently, the redistribution of the hot gas by advection tends to favour the lobe which is intersected by the majority of the captured streamlines at a given time. In the first panel of Fig. 6, for example, the hot gas is transported more efficiently into the front lobe, creating a strong front-rear asymmetry between the lobes. In the third panel, the situation is exactly the opposite and the rear lobe is more pronounced.

We point out that the continuous variations of the hot lobes and their alternating dominance are unexpected features of the heating torque which was previously thought to be strictly positive (Benítez-Llambay et al., 2015).

3.4.2 Evolving 2D and 3D streamline topology

We again explore the streamline topology and its changes related to the redistribution of the hot gas. Fig. 7 shows the midplane streamlines near the protoplanet for a selection of simulation times between $t=31.3$ and $32.15\,P_{\mathrm{orb}}$ . The former corresponds to the torque minimum, the latter to the torque maximum. The time intervals between the individual panels are not uniform but rather selected to highlight the most interesting transitions.

The sequence reveals the following features:

•

In panel (a), the streamlines are similar to the steady-state $\kappa_{\mathrm{const}}$ -disk case (bottom of Fig. 5) in several ways, mostly in the detachment of the classical horseshoe streamlines and in the existence of the retrograde streamlines captured from downstream horseshoe regions. But there are also important differences: There are none inner circulating streamlines intersecting the Hill sphere and moreover, a larger part of the captured streamlines originate in the rear horseshoe region while the front horseshoe region is almost disconnected.

•

In panel (b), the front horseshoe region becomes entirely disconnected. The captured streamlines originate exclusively in the rear horseshoe region.

•

In panel (c), also the outer circulating streamlines stop crossing the Hill sphere. Some of the streamlines that were captured in (b) now overshoot the protoplanet and make a U-turn ahead of it. The centre around which the captured streamlines enclose becomes shifted behind the protoplanet.

•

In panel (d), the front X-point moves closer to the Hill sphere and so do the front horseshoe streamlines. The rear horseshoe region becomes radially narrower and the number of U-turn streamlines overshooting the protoplanet diminishes.

•

In panel (e), the front horseshoe region reconnects with the captured streamlines.

•

In panel (f), the captured streamlines originate mostly in the front horseshoe region while the rear horseshoe region is evolving towards its disconnection, similarly to what we saw for the front horseshoe region in panels (a) and (b). The centre of the captured streamlines moves inwards from the protoplanet (and will continue to propagate ahead).

Fig. 8 shows three selected snapshots of 3D streamlines ( $t=31.3$ , $31.75$ and $32.15\,P_{\mathrm{orb}}$ ) when the oscillating torque is at its minimum (top), grows halfway towards the maximum (middle) and reaches it (bottom). Clearly, the reshaping of the streamline topology that we described for the midplane propagates in a complicated way into the vertical direction as well:

•

In the first panel, the spiraling streamlines of the vertical column are more or less centred above the protoplanet and as they rise above the midplane they penetrate the majority of the Hill sphere.

•

In the second panel, we see the overshooting rear horseshoe streamlines making their U-turns within the Hill sphere. The vertical column of captured streamlines is displaced to the rear of the Hill sphere. As the captured streamlines spiral up, the column tilts towards the Hill sphere. The outer circulating streamlines are strongly uplifted towards colatitudes above the Hill sphere.

•

In the third panel, some uplifted outer circulating streamlines penetrate into the vertical column and by this reconfiguration, the column reconnects with the front horseshoe region while the rear one starts to disconnect.

3.5 Physical processes distorting the gas flow

In previous sections, we revealed that the gas heating from an accreting protoplanet changes the topology of the flow. Perturbed streamlines have a tendency to bent towards the protoplanet and also to rise vertically. If the perturbations become strong enough, the streamlines can form a vertical spiral. In this section, we investigate the physical processes responsible for such a streamline distortion.

In Sect. 3.5.1, we theorise that the streamline distortion is a result of vorticity perturbations arising because the vigorous accretion heating renders the circumplanetary gas baroclinic. A confirmation is provided for the steady state of the $\kappa_{\mathrm{const}}$ -disk simulation in Sects. 3.5.2 and 3.5.3. Finally, Sect. 3.5.4 demonstrates that the vertical temperature gradient above the accreting protoplanet is superadiabatic and we also highlight differences between the $\kappa_{\mathrm{const}}$ -disk and $\kappa_{\mathrm{BL}}$ -disk.

3.5.1 Vorticity evolution

The spiral-like structure of the captured streamlines and the bending of nearby circulating streamlines suggests that the vorticity of the flow is modified when the protoplanet becomes hot. The vorticity can be defined via the relation

[TABLE]

where $\@vec{\omega_{\mathrm{a}}}$ is the absolute vorticity and $\@vec{\omega_{\mathrm{r}}}$ is the relative vorticity in the reference frame corotating with the protoplanet.

Evolution of $\@vec{\omega_{\mathrm{r}}}$ is described by the vorticity equation in the corotating frame (see Appendix C)

[TABLE]

where $\mathrm{D}/\mathrm{D}t$ denotes the Lagrangian derivative. In writing the equation, we neglected the effects of viscous diffusivity (large Reynolds number limit) but otherwise the equation is general.

Regarding the right-hand side terms, the first one describes the tendency of vortex tubes to become twisted due to velocity field gradients. The second one characterises the stretching or contraction of vortex tubes due to flow expansion or compression. These first two terms, usually called the twisting and stretching term, are only important if there is non-zero absolute vorticity already existing in the flow and they can cause its redistribution.

The remaining term on the right-hand side is the baroclinic term. It vanishes in barotropic flows where the pressure and density gradients are always parallel, but since our model is not barotropic, $\nabla\rho$ and $\nabla P$ can be misaligned, leading to vorticity production or destruction since their cross product can be non-zero. Because $\rho$ near the accreting protoplanet exhibits asymmetric perturbations while $P$ remains roughly spherically symmetric, one can expect non-zero baroclinic perturbations to arise. The respective non-zero vorticity then enhances circulation around a given point of the continuum, twisting the streamlines with respect to the situation unperturbed by accretion heating.

3.5.2 Baroclinic vorticity generation

It is not a priori evident which source term is the most important for the vorticity evolution in our simulations. In Fig. 9, we study the variations of $\@vec{\omega_{\mathrm{r}}}$ and the source terms of Eq. (14) along a single outer circulating streamline. We compare the situation near the cold and hot protoplanet in the $\kappa_{\mathrm{const}}$ -disk.

Downstream, before the streamline encounters the protoplanet, the situation is similar for the compared cases: the azimuthal and radial vorticity components $\omega_{\mathrm{r},\theta}$ and $\omega_{\mathrm{r},r}$ are negligible, while the polar component $\omega_{\mathrm{r},\phi}$ is non-zero and positive. The positive value of $\omega_{\mathrm{r},\phi}$ corresponds to the inherent vorticity in a flow with the Keplerian shear: Taking $v_{r}=0$ , $v_{\phi}=0$ and $v_{\theta}=\sqrt{GM/r}-\Omega r$ , one obtains $\omega_{\mathrm{r},\phi}=-0.5\Omega_{\mathrm{K}}\left(r\right)+2\Omega$ . At the same time, the source terms are zero because far from the protoplanet there are no strong velocity gradients, no compression and $\nabla\rho$ and $\nabla P$ are aligned.

To relate the vorticity variation with the streamline distortion close to the protoplanet, let us make a thought experiment, utilising the fact that the vorticity describes the tendency of the flow to circulate around some point in space. First we focus on the hot-protoplanet case which is the most important for us. Imagine an observer moving along the critical outer circulating streamline, corresponding to the outer thick black curve in the bottom right panel of Fig. 5. The observer moves with the flow, predominantly in the $-\theta$ direction. For this experiment, we dub the directions $\theta$ , $-\theta$ , $r$ , $-r$ , $\phi$ and $-\phi$ as behind, ahead, outwards, inwards, down and up, respectively. Considering only the streamlines of Fig. 5 originating at $r>r_{\mathrm{p}}$ , they initially define a plane at constant $\phi=\pi/2-0.005\,\mathrm{rad}$ and the observer propagates through that plane. When the flow reaches the Hill sphere, the observer studies the instantaneous displacement of nearby gas parcels which will correspond to the deformation of the surface defined by neighbouring streamlines.

According to bottom panels of Fig. 9, the observer measures $\omega_{\mathrm{r},\theta}>0$ when crossing the Hill sphere. Thus $\omega_{\mathrm{r},\theta}$ points against the direction of the observer’s motion, forcing the circulation in the local $(r,\phi)$ plane. The nearby gas parcels must obey the following right-hand rule: When the thumb points in the direction of $\omega_{\mathrm{r},\theta}$ , wrapping fingers determine the direction of circulation. From the observer’s point of view, an outer gas parcel will be falling downwards to the midplane and an inner gas parcel will be rising upwards. This is exactly in accordance with Fig. 5 where streamlines passing the protoplanet are uplifted from the $\phi=\pi/2-0.005\,\mathrm{rad}$ plane and the kick gets stronger with decreasing separation to the protoplanet.

$\omega_{\mathrm{r},r}$ only slightly oscillates during the Hill sphere crossing but becomes positive (although relatively small) upstream. Using again the same considerations as above, $\omega_{\mathrm{r},r}$ points outwards from the observer after crossing the Hill sphere, promoting circulation in the $(\theta,\phi)$ plane. Using the right-hand rule, a gas parcel ahead of the observer will be falling downwards and a gas parcel behind the observer will be rising upwards. In Fig. 5, this is reflected by the streamline topology when the red streamlines rising after the Hill sphere passage suddenly start to fall back towards the midplane.

As for the remaining vorticity component $\omega_{\mathrm{r},\phi}$ , it remains positive during the Hill sphere crossing, but it acquires a positive boost at first and a more prominent negative perturbation afterwards. The later diminishes the circulation related to shear in the $(r,\theta)$ plane. This is only possible if gas parcels near the observer become displaced towards trajectories with smaller shear velocities. In Fig. 5, the streamline topology indeed exhibits such a behaviour because when the red streamlines (and similarly green and yellow ones) pass the protoplanet, they are being bent towards it.

Looking at the source terms, it is obvious that the baroclinic term is responsible for perturbing $\omega_{\mathrm{r},\theta}$ and $\omega_{\mathrm{r},\phi}$ , while counteracting the twisting term contributing to $\omega_{\mathrm{r},r}$ . The importance of the baroclinic term for the flow approaching the hot protoplanet is thus confirmed. Moreover, the perturbation is indeed 3D as each of the studied components is important for the resulting streamline topology.

Comparing the hot-protoplanet case to the cold-protoplanet case, we notice that the evolution of $\omega_{\mathrm{r},\theta}$ and $\omega_{\mathrm{r},\phi}$ is roughly antisymmetric, as well as the evolution of the baroclinic source term. This is consistent with our finding that streamlines near the cold protoplanet are distorted in the opposite manner (they fall towards the midplane and deflect away from the protoplanet). In other words, the baroclinic behaviour of the protoplanet’s vicinity is reverted between the cold- and hot-protoplanet case.

3.5.3 Baroclinic region

Although we do not repeat the vorticity analysis for the remaining sets of streamlines, it is clear that features of the streamline distortion can be explained by the baroclinic generation of the vorticity. To further support this claim, we compare in Fig. 10 the map of the baroclinic term near the cold and hot protoplanet in the $\kappa_{\mathrm{const}}$ -disk. We plot the polar component of the baroclinic term $\left(\nabla\rho\times\nabla P\right)_{\phi}/\rho^{2}$ in the midplane. Since the midplane flow is effectively 2D (because $v_{\phi}=0$ in midplane), only the polar component of $\@vec{\omega_{\mathrm{r}}}$ is non-zero and thus also the polar component of the baroclinic term is the most important one.

Fig. 10 reveals that the gas is baroclinic near both the cold and hot protoplanet, but for the latter case, the map becomes approximately antisymmetric compared to the cold-protoplanet case, slightly rotated in the retrograde sense and the baroclinic region is more extended.

The existence of the baroclinic region can be explained using the isocontours of constant volume density and isobars. For a barotropic gas ( $\nabla\rho\parallel\nabla P$ ), any nearby isocontours of constant $\rho$ and $P$ should have the same shape. Wherever the isocontours depart from one another, it means that the local gradients of $\rho$ and $P$ are misaligned and that the baroclinic term is non-zero.

Looking at the cold-protoplanet case, the isobars and contours of constant density are nearly spherically symmetric. However, we notice that each density isocontour exhibits two bumps and appears to be stretched in the direction where the gas outflows from the Hill sphere. Clearly, these bumps are associated with the cold-finger perturbation that appears in the same location (see Fig. 3). Since the cold fingers are filled with overdense gas, they perturb the local density gradient, making it to point towards them. At the same time, the isobars remain approximately spherically symmetric.

When the protoplanet is hot, the cold fingers are replaced with hot underdense perturbations. Therefore, the local density gradient tends to point away from them. We can see that this is indeed true because the density isocontours do not exhibit bumps, but rather concavities across the overheated region (compare with Fig. 3). Clearly, this is the reason why the baroclinic map for the hot protoplanet looks reverted compared to the cold protoplanet.

Summarising these findings, we saw that the thermal perturbations associated either with the cold-finger effect or the heating torque make the circumplanetary region baroclinic. But the influence on the vorticity evolution is the opposite when comparing the cold- and hot-protoplanet case.

3.5.4 Vertical convection

Although the baroclinic distortion of the gas flow is a robust mechanism, it does not provide a simple explanation why the $\kappa_{\mathrm{BL}}$ -disk simulation exhibits gas instability whereas the $\kappa_{\mathrm{const}}$ -disk simulation remains stable. We now explore the vertical stability of both disks against vertical convection, considering only the hot-protoplanet limit.

To do so, we employ the Schwarzschild criterion

[TABLE]

where the subscript ‘rad’ denotes the vertical333We use the colatitude $\phi$ to study the vertical gradients because the curvature of spherical coordinates near the midplane does not significantly depart from the true vertical direction. temperature gradient found in our simulations and the subscript ‘ad’ denotes the temperature gradient that an adiabatic gas would establish. The $\nabla$ symbol stands for the logarithmic gradient $\mathrm{d}\log{T}/\mathrm{d}\log{P}$ , yielding $|\nabla_{\mathrm{ad},\phi}|=(\gamma-1)/\gamma$ for the adiabatic case.

The Schwarzschild criterion is not necessarily a universal way to determine if the disk is unstable to convection. The reason is that convective destabilisations in protoplanetary disks are opposed by diffusive effects (e.g. Held & Latter, 2018) and shear motions (Rüdiger et al., 2002). However, to our knowledge there are no convective criteria that would take into account the disk perturbation by the protoplanet, we thus use the Schwarzschild criterion for its simplicity, keeping the limitations in mind.

Fig. 11 shows the vertical velocity field and balance of the Schwarzschild criterion in the vertical direction of the $\kappa_{\mathrm{const}}$ -disk and $\kappa_{\mathrm{BL}}$ -disk. Each panel shows a different simulation time and vertical plane. For the $\kappa_{\mathrm{const}}$ -disk, we display $t=60\,P_{\mathrm{orb}}$ and the vertical plane intersecting the protoplanet’s location, whereas for the $\kappa_{\mathrm{BL}}$ -disk, each plane approximately intersects the centre around which the captured streamlines of Fig. 6 circulate at the given simulation time ( $t=31.3$ , $31.75$ , $32.2$ and $32.65\,P_{\mathrm{orb}}$ ). In other words, we choose the vertical planes where we expect the most prominent vertical outflow.

The first thing we point out is that the background differs between the studied disks. The background of the $\kappa_{\mathrm{BL}}$ -disk is slightly superadiabatic, contrary to the $\kappa_{\mathrm{const}}$ -disk. The difference arises as a result of the opacity laws. As derived by Lin & Papaloizou (1980) and Ruden & Pollack (1991), one can make a qualitative estimate for optically thick regions unperturbed by the protoplanet

[TABLE]

where $\beta$ is the power-law index of the $\kappa\propto T^{\beta}$ dependence. In the $\kappa_{\mathrm{const}}$ -disk, $\beta=0$ and thus $|\nabla_{\mathrm{rad},\phi}|/|\nabla_{\mathrm{ad},\phi}|-1\simeq-0.17$ . In the $\kappa_{\mathrm{BL}}$ -disk, the Bell & Lin (1994) opacity law in the given temperature range corresponds to water-ice grains and exhibits $\beta=2$ , therefore $|\nabla_{\mathrm{rad},\phi}|/|\nabla_{\mathrm{ad},\phi}|-1\simeq 0.67$ . These estimated values are in a good agreement with the background values of Fig. 11.

The background itself is not convective because vertical convection is usually not self-sustainable unless there is a strong heat deposition within the disk (e.g. Cabot, 1996; Stone & Balbus, 1996; Klahr et al., 1999; Lesur & Ogilvie, 2010) which, however, can be provided by the hot accreting protoplanet in our case. Indeed, looking at the $\kappa_{\mathrm{const}}$ -disk in Fig. 11 (top), the region where we previously identified the most significant temperature excess due to planetary luminosity (compare with Fig. 4) is superadiabatic and the corresponding vertical outflow can be considered convective.

In the $\kappa_{\mathrm{BL}}$ -disk simulation, the temperature gradient departs from the adiabatic one even more and additionally, the excess is no longer centred above the protoplanet itself but rather spans its vicinity. This can be seen from the varying azimuthal coordinate of the displayed vertical planes and also from the radial offset of the highly superadiabatic region in panels 3 and 5.

We summarise the section by speculating that the $\kappa_{\mathrm{BL}}$ -disk simulation becomes destabilised because the hot disturbance created by the accreting protoplanet is subject to vertical buoyant forces acting over a more extended region (compared to the $\kappa_{\mathrm{const}}$ -disk). The reason is that the hot disturbance is imposed over an already-superadiabatic background. The uplift of the material cannot be compensated for in a stationary manner and eventually, the vertical outflow becomes offset with respect to the protoplanet and starts to change its position in a cyclic manner. However, such a description is rather qualitative and precise conditions for triggering the instability should be explored in future works.

3.6 Torque oscillation vs opacity gradient

We now perform a partial exploration of the parametric space by varying the opacity gradient within the disk. The aim is twofold: First, we would like to support the claim of the previous Sect. 3.5.4 about the importance of the vertical stratification for the torque oscillations. Second, it is desirable to show that the appearance of torque oscillations in the $\kappa_{\mathrm{BL}}$ -disk is not coincidental and that it can be recovered for a wider range of parameters.

We construct six additional disk models with artificial opacity laws that i) conserve the opacity value at the protoplanet location ( $1.11\,\mathrm{cm^{2}\,g^{-1}}$ ) and ii) lead to opacity gradients which are intermediate between the $\kappa_{\mathrm{const}}$ - and $\kappa_{\mathrm{BL}}$ -disks. The latter property accounts for the highly unconstrained size distribution of solid particles in protoplanetary disks which manifests itself in a large parametric freedom of the power-law slope of the opacity profile (e.g. Piso et al., 2015).

The first additional set of disks utilises the opacity law which we dub $T$ -dependent:

[TABLE]

Similarly to the Bell & Lin (1994) opacity in the water-ice regime, it is exclusively a function of temperature. The temperature $\bar{T}$ is again azimuthally averaged to disentangle the influence of global opacity gradients (which we focus on) from those related to accretion heating of the protoplanet. We examine the values $\beta=1.5$ , $1$ and $0.5$ to span the range between $\beta=0$ ( $\kappa_{\mathrm{const}}$ -disk) and $2$ ( $\kappa_{\mathrm{BL}}$ -disk). The constant of proportionality $\kappa_{0}$ is always chosen to recover $\kappa=1.11\,\mathrm{cm^{2}\,g^{-1}}$ at $r=a_{\mathrm{p}}$ and $\phi=\pi/2$ in an equilibrium disk. We find $\kappa_{0}\simeq 0.0016$ , $0.014$ and $0.122\,\mathrm{cm^{2}\,g^{-1}}$ for the respective values of $\beta$ .

The opacity law used for the second additional set of disks, referred to as $r$ -dependent, is

[TABLE]

and we choose $\delta=3$ , $2$ and $1$ (because $\delta=4$ leads to an opacity profile similar to the $\kappa_{\mathrm{BL}}$ -disk). The motivation for choosing this purely radially dependent opacity law is to distinguish between effects caused by radial and vertical opacity gradients. The latter will not appear when $\kappa=\kappa(r)$ .

The opacity profiles of these disks in radiative equilibrium are summarised in Fig. 12 which reveals that all radial opacity gradients (top two panels) are indeed shallower compared to the $\kappa_{\mathrm{BL}}$ -disk. However, all disks with $r$ -dependent opacities have zero vertical opacity gradient by construction (bottom panel), unlike disks with $T$ -dependent opacities which vertically decrease.

The torque measurements are done the same way as in our previous experiments and the results are given in Fig. 13. For disks with $T$ -dependent opacities (top), we find that the torque oscillations appear for all investigated values of $\beta$ . The oscillation amplitude, however, linearly decreases with $\beta$ and the period becomes slightly shorter as well. In case of $r$ -dependent opacities (bottom), the torque evolution does not strongly depend on $\delta$ and exhibits marginal and vanishing oscillations, as in the $\kappa_{\mathrm{const}}$ -disk case.

Since the only qualitative difference between the disks with $T$ -dependent and $r$ -dependent opacities is in the vertical opacity gradient, our torque measurements confirm the importance of the vertical structure for the torque oscillations. The torque oscillates wherever the vertical opacity profile favours superadiabatic vertical stratification and the oscillation amplitude scales with the strength of the opacity-temperature coupling. In the absence of the vertical opacity gradient, the oscillations are not established.

3.7 Evolution of a migrating protoplanet

The previous simulations were conducted with the assumption of the fixed orbit of the protoplanet and the static torque was examined. Here we explore whether our findings can be readily applied to a dynamical case when the protoplanet is allowed to radially migrate and its semimajor axis evolves. In this section, we focus only on the $\kappa_{\mathrm{BL}}$ -disk in which we found the flow instability.

Starting from $t=30\,P_{\mathrm{orb}}$ , we release the protoplanet and run the simulation until $t=60\,P_{\mathrm{orb}}$ . Fig. 14 compares the obtained dynamical torque with the previous result of our static experiments. It is obvious that the oscillating character of the torque is retained and therefore the instability of the flow operates near a moving protoplanet as well. There are differences both in the amplitude and phase of the torque oscillations, but the mean value of the torque over the simulated period of time is $\bar{\Gamma}\simeq-1.3\times 10^{-6}\,a_{\mathrm{J}}^{2}\,\Omega_{\mathrm{J}}^{2}$ . Although this value is slightly more positive than the static heating torque ( $\bar{\Gamma}\simeq-6.3\times 10^{-6}\,a_{\mathrm{J}}^{2}\,\Omega_{\mathrm{J}}^{2}$ ), it is almost the same as the torque acting on the cold protoplanet ( $\bar{\Gamma}\simeq-1.6\times 10^{-6}\,a_{\mathrm{J}}^{2}\,\Omega_{\mathrm{J}}^{2}$ ). We thus confirm that in the $\kappa_{\mathrm{BL}}$ -disk, the heating torque does not add any considerable positive contribution to the mean torque but makes it strongly oscillating instead.

Bottom panel of Fig. 14 shows the actual evolution of the protoplanet’s semimajor axis. On average, the protoplanet slowly migrates inwards but this drift is not smooth. The protoplanet exhibits fast periodic inward and outward excursions on an orbital time scale. The migration rate of these individual excursions (not to be confused with the mean migration rate stated above) is $\dot{a}\sim(10^{-3}\,\mathrm{au})/P_{\mathrm{orb}}$ .

Finally, we note that the mean migration rate is not constant, which also corresponds to the varying offset of the dynamical torque with respect to the static torque in Fig. 14. Probably, the unstable gas distribution around the protoplanet is further affected by the protoplanet’s radial drift.

4 Discussion

This section discusses the applicability of our results, draws links to some previous studies and speculates about possible implications for planet formation. Additionally, we outline possibilities for future work. Although most of them are beyond the scope of this paper, we at least present several additional simulations in Appendix B to verify the correctness of our model (we study the impact of the luminosity increase, grid resolution, opacity treatment and computational algorithm on the evolution of torque oscillations).

4.1 Relation to previous works

Our results revealed unexpected perturbations of the gas flow in the vicinity of the protoplanet undergoing strong accretion heating, one of them being the vertical outflow. The resolution of our simulations was originally tailored for studying the torque and may lack some accuracy close to the protoplanet where the outflow occurs. Follow-up studies should thus utilise simulations with high resolution of the protoplanetary envelope, similarly to e.g. Tanigawa et al. (2012); Fung et al. (2015); Ormel et al. (2015); Lambrechts & Lega (2017).

It is likely that a critical combination of the opacity $\kappa$ and luminosity $L$ for a given planetary mass $M_{\mathrm{p}}$ exists for which the vertical outflow is triggered, similarly to the dependencies of the heating torque explored by Benítez-Llambay et al. (2015). Popovas et al. (2018, 2019) studied the stability of the circumplanetary envelope during pebble accretion and found that the gas within the Bondi sphere exhibits 3D convective motions, assuming444Greater number of parameters were discussed in Popovas et al. (2018, 2019) but we quote those closest to this paper. $L\simeq 1.4\times 10^{26}\,\mathrm{erg}\,\mathrm{s}^{-1}$ , $M_{\mathrm{p}}=0.95\,M_{\oplus}$ and $\kappa=1\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$ . Lambrechts & Lega (2017) also explored a set of ( $L$ , $\kappa$ , $M_{\mathrm{p}}$ ) parameters and its impact on the structure of the circumplanetary envelope. They found the inner region of the envelope to depart from hydrostatic equilibrium when the luminosity exceeds $L=10^{27}\,\mathrm{erg}\,\mathrm{s}^{-1}$ around a $M_{\mathrm{p}}=5\,M_{\oplus}$ core within a $\kappa=1\,\mathrm{cm}^{2}\,\mathrm{g}^{-1}$ environment, but they did not identify any vertical outflow. The outflow in our simulations appears for higher $L$ and smaller $M_{\mathrm{p}}$ compared to Lambrechts & Lega (2017), we can therefore assume that our parameters cross the critical ones.

Regarding the baroclinic perturbations, although they are known to produce vortical instabilities in protoplanetary disks (Klahr & Bodenheimer, 2003; Petersen et al., 2007a, b; Lesur & Papaloizou, 2010; Raettig et al., 2013; Barge et al., 2016), they have been rarely considered in relation to hot protoplanets. For example, Owen & Kollmeier (2017) claim that that hot protoplanets can excite large-scale baroclinic vortices but we do not identify any of those in our simulations. Instead, we find baroclinic perturbations to be responsible for 3D distortion of the gas flow near the protoplanet.

4.2 Implications for the formation of planetary systems

Although the heating torque has been previously thought to be strictly positive and also efficient in high-opacity locations of protoplanetary disks, our paper shows that it can exhibit more complicated behaviour if the temperature dependence of the disk opacity is taken into account. We identified an oscillatory mode of the heating torque in the disk region with $\kappa\propto T^{2}$ and we demonstrated that it can operate even for shallower dependencies such as $\kappa\propto T^{0.5}$ , albeit with a decreased amplitude. This behaviour resembles the nature of baroclinic and convective disk instabilities which usually operate in the most opaque regions with the steepest entropy gradients but become less effective elsewhere (e.g. Pfeil & Klahr, 2019).

It is worth noting that our simulations neglected the effect of stellar irradiation. Stellar-irradiated disks tend to have vertical temperature profiles that are closer to being isothermal (e.g. Flock et al., 2013), unlike disk models used in this work which have rather adiabatic or slightly superadiabatic vertical temperature gradients. On the other hand, even the irradiated disks often contain shadowed regions protected against stellar irradiation, for example behind the puffed-up inner rim (e.g. Dullemond et al., 2001) or between the viscously heated inner disk and a flared irradiated outer disk (e.g. Bitsch et al., 2013). In such regions, the oscillatory migration could still operate.

For the aforementioned reasons, it is likely that transition zones might exist in protoplanetary disks, separating regions where protoplanets migrate under the influence of the standard positive heating torque and where they undergo oscillatory migration. Migration at the edges of such zones could be convergent, leading to a pileup of protoplanets.

5 Conclusions

By means of 3D radiation-hydrodynamic simulations, we investigated the heating torque (Benítez-Llambay et al., 2015) acting on a luminous $3\,M_{\oplus}$ protoplanet heated by accretion of solids. The aim was to compare the torque evolution and physics in a disk with non-uniform opacities (Bell & Lin, 1994) with the outcome of a constant-opacity simulation.

We discovered that the gas flow near the protoplanet is perturbed by two mechanisms:

The gas advected past the protoplanet becomes hot and underdense. Consequently, a misalignment is created between the gradients of density and pressure within the protoplanet’s Hill sphere. The baroclinic term of the vorticity equation ( $\sim$$\nabla\rho\times\nabla P$ ) then becomes non-zero and modifies the vorticity of the flow. 2. 2.

The efficient heat deposition in the midplane makes the vertical temperature gradient superadiabatic, thus positively enhancing vertical gas displacements.

The streamline topology exhibits a complex 3D distortion. The most important feature are spiral-like streamlines rising vertically above the hot protoplanet, forming an outflow column of gas escaping the Hill sphere.

In the constant-opacity disk, the vertical outflow is centralised above the protoplanet, it temporarily captures streamlines from both horseshoe regions and such a state is found to be stationary over the simulation time scale. The distribution of the hot gas then remains in accordance with findings of Benítez-Llambay et al. (2015), having a two-lobed structure, and so does the resulting positive heating torque.

In the disk with non-uniform opacity $\kappa\propto T^{2}$ (typically outside the water-ice line), we find the superadiabatic temperature gradient to be steeper and the distorted gas flow to be unstable. The vertical spiral flow becomes offset with respect to the protoplanet and periodically changes its position, spanning the edge of the Hill sphere in a retrograde fashion. Its motion is followed by the underdense gas and the resulting heating torque strongly oscillates in time. The interplay can be characterised by the following sequence:

A stage when most of the captured streamlines originate in the rear horseshoe region and their spiral-like structure is offset ahead of the protoplanet. Therefore the hot gas cumulates ahead of the protoplanet, a dominant underdense lobe is formed there and the torque becomes negative, reaching the minimum of its oscillation. 2. 2.

A stage when the front horseshoe region becomes completely isolated from the captured streamlines. Some of the rear horseshoe streamlines start to overshoot the protoplanet and make U-turns ahead of it. At the same time, the spiral-like structure recedes behind the protoplanet. The lobe from stage 1 starts to decay while a rear lobe starts to grow and the torque changes from negative to positive. 3. 3.

Antisymmetric situation to stage 1, when most of the captured streamlines originate in the front horseshoe region, the dominant lobe trails the protoplanet and the torque is positive, reaching the maximum of its oscillation. 4. 4.

Antisymmetric situation to stage 2, when the torque decreases from positive to negative and the cycle repeats.

Such an advective redistribution of the hot underdense gas is sustained over the simulation time scale.

We also studied the dependence of the torque oscillations on the opacity gradient and found that they can appear even for $\kappa\propto T^{0.5}$ , although their amplitude linearly decreases with the power-law slope of the $\kappa(T)$ dependence. We also demonstrated that the oscillations would vanish in a disk with zero vertical opacity gradient.

If the protoplanet is allowed to migrate, its mean migration rate is nearly unaffected but the radial drift is not smooth. It is rather oscillatory, consisting of brief inward and outward radial excursions with a characteristic rate $\dot{a}\sim 10^{-4}\,\mathrm{au}\,\mathrm{yr}^{-1}$ .

We discussed possible implications of the oscillating heating torque for planet formation and pointed out that it can affect the global evolution of hot migrating low-mass protoplanets. During their migration through disk regions with varying opacities it might be possible for the protoplanets to switch between the standard positive heating torque of Benítez-Llambay et al. (2015) and the positive/negative torque oscillations discovered in this paper.

Acknowledgements.

The work of OC has been supported by Charles University (research program no. UNCE/SCI/023; project GA UK no. 624119; project SVV-260441), by the Grant Agency of the Czech Republic (grant No. 18-06083S) and by The Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations project IT4Innovations National Supercomputing Centre LM2015070. Access to computing and storage facilities owned by parties and projects contributing to the National Grid Infrastructure MetaCentrum, provided under the programme “Projects of Large Research, Development, and Innovations Infrastructures” (CESNET LM2015042), is greatly appreciated. OC would like to thank M. Brož and B. Bitsch for useful discussions and also to E. Lega who kindly provided her numerical code fargoca for comparison simulations. The authors are grateful to an anonymous referee whose insightful ideas helped to improve the manuscript.

Appendix A Modifications of the numerical scheme

We implemented Eqs. (3) and (4) into fargo3d using their discrete form derived by Bitsch et al. (2013) (see their Appendix B). We introduce a minor modification of the numerical scheme, which allows to account for all the source terms in a single substep. Using the same notation as in Bitsch et al. (2013), the relation between the temperature and radiative energy at $t+\Delta t$ is

[TABLE]

and we redefine

[TABLE]

There are two changes with respect to Bitsch et al. (2013). First, the compressional heating is included via the last term in the denominator of Eqs. (20) and (21). Second, the $Q_{\epsilon\mathrm{-indep}}$ term is a sum of all heat sources that do not depend on $\epsilon$ (or $T$ ); in our case $Q_{\epsilon\mathrm{-indep}}=Q_{\mathrm{visc}}+Q_{\mathrm{art}}+Q_{\mathrm{acc}}$ . We note that when necessary, the stellar irradiation term can also be easily included in $Q_{\epsilon\mathrm{-indep}}$ but it is neglected in this work.

Appendix B Supporting simulations

In this appendix, we summarize several additional simulations designed to confirm the robustness of our conclusions.

B.1 Simulation with a smoothly increasing luminosity of the protoplanet

In our main simulations, we usually start the phase with accretion heating abruptly, by instantaneously increasing the luminosity of the protoplanet from $L=0$ to the value corresponding to the mass doubling time $\tau=100\,\mathrm{kyr}$ . Such a sudden appearance of a strong heat source might produce an unexpected behaviour and instabilities by itself. In order to exclude any undesirable behaviour, we repeat the $\kappa_{\mathrm{BL}}$ -disk simulation with accretion heating of the protoplanet, but now we linearly increase $L$ from zero at $t=30\,P_{\mathrm{orb}}$ to its maximal value over the time interval of $10\,P_{\mathrm{orb}}$ .

Fig. 15 shows the measured torque evolution. Clearly, the gradual increase of $L$ has no impact on the final character of the torque, the oscillations related to the flow reconfigurations inevitably appear.

B.2 Simulation with an increased azimuthal resolution

The resolution in our simulations is motivated by works of Lega et al. (2014) and Eklund & Masset (2017). Although we do not perform extended convergence tests of our own, the resolution should be sufficient to recover a realistic torque value and also a realistic advection-diffusion redistribution of the hot gas near the protoplanet.

However, once the circumplanetary flow becomes unstable, it is no longer clear if the chosen resolution is sufficient. For example, one might argue that the coverage of the Hill sphere by the grid cells in the azimuthal direction is too poor. Here we present an experiment in which we double the number of the grid cells in the azimuthal direction in order to have the same coverage of the Hill sphere in all directions. We perform the $\kappa_{\mathrm{BL}}$ -disk simulation again, however, we shorten the phase without accretion heating to $10\,P_{\mathrm{orb}}$ and the phase with accretion heating to $5\,P_{\mathrm{orb}}$ .

The result is shown in Fig. 16 and demonstrates that the toque oscillations are recovered even when the increased resolution is used. However, the amplitude of the torque oscillations slightly changes, implying that the resolution dependence should be explored more carefully in the future.

B.3 Simulation with the unmodified Bell & Lin opacity table

The simulations of the $\kappa_{\mathrm{BL}}$ -disk presented in this paper are performed with a simplified opacity law of Bell & Lin (1994) (explained in Sect. 2.2). Here we test how the results change if the unmodified opacity law $\kappa_{\mathrm{BL}}^{\mathrm{full}}$ is used and the dependence on the local values of $T$ and $\rho$ is retained.

Fig. 17 compares the torque evolution in our standard $\kappa_{\mathrm{BL}}$ -disk with that in a $\kappa_{\mathrm{BL}}^{\mathrm{full}}$ -disk. Clearly, the instability occurs in both disks, regardless of whether or not the input values for the opacity function are azimuthally averaged. The only difference is in the torque amplitude which is larger in the $\kappa_{\mathrm{BL}}^{\mathrm{full}}$ -disk.

The increased amplitude occurs because if $T$ locally rises, so does the material opacity ( $\kappa_{\mathrm{BL}}^{\mathrm{full}}\propto T^{2}$ in the given disk region). Subsequently, the radiative cooling of the hot gas becomes less efficient and $T$ rises even more. As a result, the underdense perturbations related to any temperature excess are more pronounced, leading to the larger amplitude of the torque oscillations.

B.4 Code comparison

With the aim to confirm that our implementation of the energy equations in fargo3d is correct and also that the instability of the circumplanetary flow does not arise due to numerical artefacts, we present a comparison simulation obtained with an independent and well-tested radiation hydrodynamic code fargoca (Lega et al. 2014). The simulation is performed with the unmodified $\kappa_{\mathrm{BL}}^{\mathrm{full}}$ opacity law of Bell & Lin (1994).

Fig. 18 compares the torque evolution found using our code with the one obtained with fargoca. We can see that the converged torque for the cold protoplanet is in a satisfactory agreement.

When the accretion heating is initiated, the same oscillatory trend is observed with both codes. The curves overlap at first; later they start to depart in terms of the oscillation phase. However, the amplitude remains the same.

Since the converged torques are in agreement and the instability is recovered, we conclude that the differences that we identify in Fig. 18 arise only because our radiation module in fargo3d relies on a slightly different numerical scheme (see Appendix A) compared to fargoca.

Appendix C Vorticity equation in the corotating frame

For reader’s convenience, we provide a step-by-step derivation of the vorticity equation. Starting with Eq. (2), we apply the curl on both sides. In the following, we neglect the viscous term (large Reynolds number limit).

The following identities of the vector calculus will be utilised:

[TABLE]

where the last identity holds for scalar functions that are at least twice continuously differentiable.

The curl of the advection term yields

[TABLE]

where we used Eq. (22) in writing the first equality, Eq. (25) to remove the $\sim$$\@vec{v}\cdot\@vec{v}$ term and we defined the relative vorticity in the corotating frame $\@vec{\omega_{\mathrm{r}}}=\nabla\times\@vec{v}$ . Using Eqs. (23) and (24), we further obtain

[TABLE]

The curl of the pressure term leads to

[TABLE]

because $\nabla\times\left(\nabla P\right)=\@vec{0}$ (Eq. 25).

When dealing with the gravitational term, it is useful to realise that the centrifugal acceleration can be expressed as $\@vec{\Omega}\times\left(\@vec{\Omega}\times\@vec{r}\right)=\nabla\Phi_{\mathrm{c}}$ , with the centrifugal potential $\Phi_{\mathrm{c}}=-\frac{1}{2}r_{\perp}^{2}\Omega^{2}$ . The curl of a combined force term, $\nabla\times\left[\nabla\left(\Phi+\Phi_{\mathrm{c}}\right)\right]$ , is zero owing to Eq. (25).

Finally, we take the curl of the Coriolis acceleration:

[TABLE]

where we removed terms $\sim$$\nabla\cdot\@vec{\Omega}$ , $\sim$$\nabla\@vec{\Omega}$ because $\@vec{\Omega}$ in our simulations is constant.

Recollecting all the terms, we can write the relative vorticity equation

[TABLE]

where we defined the absolute vorticity in the inertial frame $\@vec{\omega_{\mathrm{a}}}=\@vec{\omega_{\mathrm{r}}}+2\@vec{\Omega}$ .

Bibliography83

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Barge et al. (2016) Barge, P., Richard, S., & Le Dizès, S. 2016, A&A, 592, A 136
2Baruteau et al. (2011) Baruteau, C., Fromang, S., Nelson, R. P., & Masset, F. 2011, A&A, 533, A 84
3Baruteau & Masset (2008) Baruteau, C. & Masset, F. 2008, Ap J, 672, 1054
4Bell & Lin (1994) Bell, K. R. & Lin, D. N. C. 1994, Ap J, 427, 987
5Benítez-Llambay et al. (2015) Benítez-Llambay, P., Masset, F., Koenigsberger, G., & Szulágyi, J. 2015, Nature, 520, 63
6Benítez-Llambay & Masset (2016) Benítez-Llambay, P. & Masset, F. S. 2016, Ap JS, 223, 11
7Bitsch et al. (2013) Bitsch, B., Crida, A., Morbidelli, A., Kley, W., & Dobbs-Dixon, I. 2013, A&A, 549, A 124
8Brasser et al. (2018) Brasser, R., Matsumura, S., Muto, T., & Ida, S. 2018, Ap J, 864, L 8