Oscillatory migration of accreting protoplanets driven by a 3D distortion of the gas flow
AAstronomy & Astrophysics manuscript no. paper c (cid:13)
ESO 2019April 30, 2019
Oscillatory migration of accreting protoplanetsdriven by a 3D distortion of the gas flow
O. Chrenko and M. Lambrechts Institute of Astronomy, Charles University in Prague, V Holešoviˇckách 2, CZ–18000 Prague 8, Czech Republice-mail: [email protected] Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University, Box 43, 22100 Lund, SwedenReceived ??; accepted ??
ABSTRACT
Context.
The dynamics of a low-mass protoplanet accreting solids is influenced by the heating torque, which was found to suppressinward migration in protoplanetary disks with constant opacities.
Aims.
We investigate the di ff erences 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 thelocal density and pressure gradients. As a result, the 3D gas flow is perturbed and some of the streamlines form a retrograde spiralrising above the protoplanet. In the constant-opacity disk, the perturbed flow reaches a steady state and the underdense gas responsiblefor the heating torque remains distributed in accordance with previous studies. If the opacity is non-uniform, however, the di ff erencesin the disk structure can lead to more vigorous streamline distortion and eventually to a flow instability. The underdense gas developsa one-sided asymmetry which circulates around the protoplanet in a retrograde fashion. The heating torque thus strongly oscillates intime 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 alternatingrapid inward and outward migration. We speculate that transitions between the positive and oscillatory heating torque may occur inspecific 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 disksis a key evolutionary step in formation of each planetary sys-tem. Low-mass protoplanets, incapable of gap opening (Cridaet al. 2006), undergo Type I migration under the influence ofthe gravitational torques exerted by the Lindblad spiral wakes(Goldreich & Tremaine 1979; Ward 1986) and by the gasin their corotation region (Ward 1991; Masset 2002; Tanakaet al. 2002; Paardekooper & Mellema 2006; Baruteau & Masset2008; Paardekooper & Mellema 2008; Masset & Casoli 2009;Paardekooper & Papaloizou 2009; Baruteau et al. 2011). TypeI migration depends in a complicated way on the disk structureand thermophysics (e.g. Kley & Crida 2008; Paardekooper et al.2010, 2011; Lega et al. 2015). Its detailed understanding is there-fore essential for the creation of realistic population synthesismodels (e.g. Coleman & Nelson 2016).It has been recently discovered that low-mass protoplanetsevolving in radiative disks are subject to thermal torques (Legaet al. 2014; Benítez-Llambay et al. 2015; Masset & VelascoRomero 2017; Masset 2017) related to the thermal perturbationsinduced by a protoplanet in its vicinity. If the protoplanet itselfis cold and non-luminous (Lega et al. 2014), the gas which ar-rives into its potential well becomes heated mostly as a result ofcompression (by means of the thermodynamic ‘ P d V ’ term). Thearising temperature excess becomes smoothed out by the radia-tive transfer, so when the gas leaves the high pressure region it lacks some of its internal energy – compared to the state beforethe compression – and therefore becomes colder and overdense.Two overdense lobes appear along the streamlines outflowingfrom the Hill sphere and their asymmetry makes the total torquefelt by the protoplanet more negative, enhancing the inward mi-gration. The process is known as the cold-finger e ff ect (Legaet al. 2014).In the opposite limit, the protoplanet is hot, as a result of thesolid material deposition during its formation (e.g. by pebble ac-cretion; Ormel & Klahr 2010; Lambrechts & Johansen 2012). Insuch a case, the luminous protoplanet acts as a local heat sourcefor the surrounding gas. Once the gas is heated, it becomes un-derdense compared to the situation without accretion heating.Benítez-Llambay et al. (2015) performed 3D radiation hydro-dynamic simulations with the assumption of the constant diskopacity and found that the hot protoplanet on a fixed circular or-bit creates two underdense lobes of gas, again associated withthe outflow from the Hill sphere. The rear lobe (positioned be-hind the protoplanet outwards from its orbit) is dominant, there-fore there is an overabundance of gas ahead of the protoplanetand the resulting torque becomes more positive, supporting out-ward migration. The positive enhancement was named the heat-ing torque. It was proposed to be an additional mechanism (alongwith 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 & Chiang2017; Brasser et al. 2018; Miranda & Lai 2018; McNally et al. Article number, page 1 of 19 a r X i v : . [ a s t r o - ph . E P ] A p r & A proofs: manuscript no. paper ffi cient inward migration (Ko-rycansky & Pollack 1993; Ward 1997; Tanaka et al. 2002).Moreover, the heating torque has important dynamical con-sequences for migrating protoplanets (Brož et al. 2018; Chrenkoet al. 2018) because it can excite orbital eccentricities and incli-nations by means of the hot-trail e ff ect (Eklund & Masset 2017;Chrenko et al. 2017), which counteracts the otherwise e ffi cienteccentricity and inclination damping by waves (Tanaka & Ward2004; Cresswell et al. 2007).Nevertheless, the heating torque has not been extensivelystudied in 3D radiative disks with non-uniform opacities. How-ever, realistic opacity functions (Bell & Lin 1994; Semenov et al.2003; Zhu et al. 2012) are of a great significance for the diskstructure and planet migration (e.g. Kretke & Lin 2012; Bitschet al. 2013) and this work will show that the heating torque isa ff ected as well.In this paper, we reinvestigate the thermal torques acting ona low-mass protoplanet, with a special emphasis on the heat-ing torque in a disk with non-uniform opacity. We examine thestreamlines near the protoplanet and point out the importance oftheir 3D distortion for redistribution of the hot underdense gasresponsible for the heating torque.
2. Model
We consider a protoplanetary system consisting of a central pro-tostar surrounded by a disk of coupled gas and dust in which asingle protoplanet is embedded. The fluid part of the disk modelaccounts only for the gas, assuming the dust is a passive tracerthat acts as the main contributor to the material opacity.
The disk is described using Eulerian hydrodynamics on a stag-gered spherical mesh centred on the protostar. The spherical co-ordinates consist of the radial distance r , azimuthal angle θ andcolatitude φ . Our model is built on top of the hydrodynamicmodule of fargo d (Benítez-Llambay & Masset 2016), whichsolves the equations of continuity and momentum of a fluid ∂ρ∂ t + ( v · ∇ ) ρ = − ρ ∇ · v , (1) ∂ v ∂ t + ( v · ∇ ) v = − ∇ P ρ −∇ Φ+ ∇ · T ρ − [2 Ω × v + Ω × ( Ω × r )] , (2)where the individual quantities stand for the volume density ρ ,time t , flow velocity vector v (with the radial, azimuthal and ver-tical components v r , v θ and v φ ), pressure P , gravitational poten-tial Φ of the protostar and the protoplanet, viscous stress tensor T (Benítez-Llambay & Masset 2016), radius vector r and angu-lar velocity vector Ω which is non-zero because we work in thereference frame corotating with the protoplanet.To account for the energy production, dissipation and trans-port in the disk, we implement the two-temperature energy equa-tions for the gas and the radiation field according to the formu-lation of Bitsch et al. (2013): ∂ E R ∂ t + ∇ · F = ρκ P (cid:104) σ T − cE R (cid:105) , (3) Public version of the code is available at http://fargo.in2p3.fr/ . ∂(cid:15)∂ t + ( v · ∇ ) (cid:15) = − P ∇ · v − ρκ P (cid:104) σ T − cE R (cid:105) + Q visc + Q art + Q acc , (4)where E R is the radiative energy, F the radiation flux, κ P thePlanck opacity, σ the Stefan-Boltzmann constant, T the gas tem-perature, c the speed of light, (cid:15) the internal energy of the gas, Q visc the viscous heating term (Mihalas & Weibel Mihalas 1984), Q art describes the heating due to the artificial viscosity (Stone &Norman 1992) and Q acc the heat released when the protoplanetis accreting. Stellar irradiation is neglected in this paper for sim-plicity although our code is capable of including it as well.The state equation of the ideal gas is used: P = ( γ − (cid:15) = ( γ − ρ c V T , (5)where γ is the adiabatic index and c V is the specific heat at con-stant volume, which can be expressed as c V = R / ( µ ( γ − R is the universal gas constant and µ is the mean molecu-lar weight.The flux-limited di ff usion approximation (FLD; Levermore& Pomraning 1981) is adopted to obtain a closure relation forthe radiation flux F = − λ lim c ρκ R ∇ E R , (6)where κ R is the Rosseland opacity and λ lim is the flux limiterof Kley (1989). For the opacities, we assume that the Planckand Rosseland means are similar enough to be replaced with asingle value κ . This is a valid assumption in the cold regions ofprotoplanetary 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 L = GM p ˙ M p R p = GM R p τ , (7)where G is the gravitational constant, M p is the mass of the pro-toplanet, ˙ M p is its accretion rate, and R p is the protoplanet’s ra-dius. In writing the second equality, we introduce the mass dou-bling time τ = M p / ˙ M p , which is a free parameter that controlsthe accretion rate in our model (see Sect. 2.2).We assume that the radiation flux from the protoplanet iscompletely absorbed by the optically thick gas in the eight gridcells enclosing the protoplanet (Benítez-Llambay et al. 2015;Eklund & Masset 2017; Lambrechts & Lega 2017). The accre-tion heating, which is non-zero only within these cells, is thensimply Q acc = LV , (8)where V is the total volume of the heated cells.The disk evolves in the combined gravitational potential ofthe central protostar and the embedded protoplanet: Φ = − GM (cid:63) r − GM p d f sm , (9)where M (cid:63) is the mass of the protostar and d is the cell-protoplanet distance. The planetary potential is smoothed toavoid numerical divergence at the protoplanet location ( d = Article number, page 2 of 19. Chrenko and M. Lambrechts: Oscillatory migration of accreting protoplanets driven by a 3D distortion of the gas flow using the tapering cubic-spline function f sm of Klahr & Kley(2006): f sm = , ( d > r sm ) , (cid:32) dr sm (cid:33) − (cid:32) dr sm (cid:33) + dr sm , ( d ≤ r sm ) (10)where the smoothing length r sm is a fraction of the protoplanet’sHill sphere radius (see Sect. 2.2). Orbital evolution of the pro-toplanet is tracked using the ias
15 integrator (Rein & Spiegel2015) from the rebound package (Rein & Liu 2012), which weinterfaced with fargo d .Although fargo d is an explicit hydrodynamic solver, imple-menting Eqs. (3) and (4) in an explicit form would lead to a veryrestrictive Courant condition on the largest allowed time steplength. To avoid such a time step limitation, it is advantageousto solve the energy equations in an implicit form. We thus fol-low 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 r ≡ A x − b , where A is the ma-trix of the linear system, x is the solution vector and b is theright-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 (cid:107) r (cid:107) / (cid:107) b (cid:107) < − , where the norms are calculated in the L space. Initial conditions describing the disk are adopted from Kley et al.(2009) and Lega et al. (2014). The surface density follows thepower-law function
Σ = r / (1 au)) − . g cm − . The initial Σ isconverted into the initial ρ assuming the disk is vertically isother-mal. The velocity field is given by balancing the gravitationalacceleration from the protostar, the acceleration due to the pres-sure gradient and the centrifugal acceleration. We assign con-stant kinematic viscosity ν = cm s − to the gas to mimicthe angular momentum transport in the disk driven by physicale ff ects outside the scope of this study (e.g. the turbulent eddyviscosity). The initial aspect ratio is h = H / r = .
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 γ = .
43 and the mean molecular weightis µ = . − , corresponding to the solar mixture of the hy-drogen and helium.We will assume a single embedded super-Earth with the mass M p = M ⊕ , orbiting on a circular orbit at the Jupiter’s he-liocentric distance a p ≡ a J = . Ω p ≡ Ω J = (cid:113) GM (cid:63) / a and the smoothing lengthfor the tapering function of the gravitational potential f sm is r sm = . R H , where R H = a p ( M p / (3 M (cid:63) )) / is the Hill sphereradius. We will study both cases without and with planetary ac-cretion, corresponding to the cold- and hot-protoplanet limit, re-spectively. In the simulations with the accretion, we assume themass doubling time τ =
100 kyr, which is a value within the Public version of the code is available at https://rebound.readthedocs.io/en/latest/ . opa c i t y κ [ c m g - ] protoplanetlocation κ const κ BL a s pe c t h = H / r
60 90 120 150 t e m pe r a t u r e T [ K ] en t r op y P / ρ γ [ c . u .] radial distance r [au] Fig. 1.
Radial profiles of the characteristic quantities in the κ const -disk(dashed blue curves) and the κ BL -disk (solid black curves) after theirnumerical relaxation over t = P orb . From top to bottom, the panelsshow the opacity κ , aspect ratio h = H / r , temperature T and entropy S = P /ρ γ . Midplane values are displayed (or used to derive h ). The dottedvertical line indicates the protoplanet’s orbital distance a p = . range of the expected pebble accretion rates (Lambrechts & Jo-hansen 2014; Chrenko et al. 2017). The resulting luminosity ofthe protoplanet is L (cid:39) . × erg s − . Using the initial disk parameters, we setup two fiducial diskmodels that only di ff er in the material opacity function. Theopacity is either constant or follows a slightly modified (see be-low) prescription of Bell & Lin (1994). To distinguish betweenthe models, we use the abbreviations κ const -disk and κ BL -disk,respectively. The exact value of the opacity in the κ const -disk istuned to be the same as in the κ BL -disk at the location of theprotoplanet a p , leading to κ const ≡ κ BL ( r = a p ) = .
11 cm g − .The unmodified opacity from Bell & Lin (1994), which wedenote κ fullBL ( ρ ( r , θ, φ ) , T ( r , θ, φ )), is a fitting law that sets κ asa function of the local gas density and temperature. The tablespans several regimes corresponding to the presence (or absence)of dust or molecular species dominant in protoplanetary disks.Using κ fullBL in our simulations with the accretion heating of theprotoplanet could cause strong local opacity gradients, becausewe expect the temperature perturbations to reach ∼ K at dis-tances of one cell size from the centre of the protoplanet.In practice, we apply the Bell & Lin’s opacity law in a sim-plified way, using κ BL ( ¯ ρ ( r , φ ) , ¯ T ( r , φ )), where the bared quantitiesare azimuthally computed arithmetic means and the dependenceon the θ -coordinate is therefore dropped. This helps us to dis-tinguish the e ff ects caused by the global structure of the diskfrom those related to the local κ - T - ρ feedback. To justify our Article number, page 3 of 19 & A proofs: manuscript no. paper approach, we verify in Appendix B that the unmodified opacitytable κ fullBL of Bell & Lin (1994) does not change our conclusions.The fiducial κ const - and κ BL -disks are discussed throughoutthe majority of the paper, with the exception of Sect. 3.6 wherewe study the dependence of the heating torque on the opacitygradient within the disk. Migration of low-mass protoplanets critically depends on thegrid resolution, we therefore combine the well-established diskextent and resolution from Lega et al. (2014) (in the azimuthaland vertical direction) and Eklund & Masset (2017) (in the ra-dial direction). The disk radially stretches from r min = .
12 au to r max = .
28 au and is resolved by 512 rings. We prevent any ver-tical motions of the protoplanet in our simulations by assumingthat the solution is symmetric with respect to the midplane. Oneof the disk boundaries in the colatitude is therefore located at themidplane ( φ = π/ ◦ . The colatitude is resolved by 64 zones. In the azimuth, weuse only 1 sector in our preparatory simulations without the pro-toplanet and 1382 sectors in our simulations with the embeddedprotoplanet. The resulting local resolution is 8 cells per R H in the r - and φ -directions and 3 cells per R H in the θ -direction (see alsoAppendix B.2 for a simple resolution test).The azimuthal boundary conditions are periodic for all prim-itive quantities. The radial boundary conditions are symmetricfor ρ , (cid:15) , and v φ and reflecting for v r . E R is set to a zero gra-dient and v θ is extrapolated using the same radial dependenceas for the Keplerian rotation velocity. The boundary conditionsin colatitude are symmetric for ρ , (cid:15) , v r and v θ , and reflectingfor v φ . E R is symmetrised at the midplane boundary and set to a R T at the remaining boundary in the colatitude, where a R isthe radiation constant and T ≡
3. Simulations
Since we use a non-isothermal equation of state and we alsoaccount for the energy production and transfer, the parametricsetup of the disk, which was discussed so far, is not station-ary. Therefore, before the simulations with the embedded pro-toplanet are conducted, we let both disks ( κ const -disk and κ BL -disk) numerically evolve over the time scale 100 P orb , where P orb = π/ Ω J .The equilibrium state after the relaxation is characterized byFig. 1. We plot the radial profiles of the midplane temperature T ,opacity κ , entropy S = P /ρ γ and aspect ratio h = H / r , where thepressure scale height is determined as H = c s / Ω K and the soundspeed as c s = (cid:112) γ P /ρ .The constant opacity of the κ const -disk makes all the remain-ing radial dependencies rather shallow. The aspect ratio is al-most flat, only slightly increasing with the radial distance. In the κ BL -disk, on the other hand, the opacity has a peak near 3 . ffi ciency of the disk cooling in-creases and simultaneously, the e ffi ciency of the viscous heatingdiminishes with the dropping relative velocity of the shearinglayers. Consequently, the aspect ratio radially decreases as theenergy budget is not su ffi cient to pu ff up the disk. In such a disk, -5e-05-4e-05-3e-05-2e-05-1e-05 0 1e-05 2e-05 3e-05 4e-05 0 10 20 30 40 50 60 t o t a l t o r que Γ [ a J Ω J ] no. of orbits [ P orb ]Hot protoplanetCold protoplanet κ const κ BL Fig. 2.
Temporal evolution of the total torque exerted on the protoplanet( M p = M ⊕ ) by the κ const -disk (dashed blue curve) and the κ BL -disk(solid black curve). For t ≤ P orb , the protoplanet is not accreting. Itsaccretion and the respective heating are activated for t > P orb , as alsohighlighted by a yellow background. While the blue curve evolves asexpected, i.e. gains a positive boost when the accretion heating is initi-ated (Benítez-Llambay et al. 2015), the black curve does not convergeand exhibits strong oscillations between positive and negative valuesinstead. The red horizontal dotted line shows the mean value of the os-cillating black curve and demonstrates that the heating torque makes thetotal torque more negative in this case. T and S profiles radially decrease more steeply compared to the κ const -disk. Starting from the relaxed state of the disks, we copy the hydrody-namic quantities in the azimuthal direction to expand the resolu-tion from a single sector to the desired 1382 sectors. The proto-planet is inserted and we simulate 30 P orb of evolution while ne-glecting any accretion and accretion heating of the protoplanet.The aim of this part of the simulation is to let the disk adjust tothe presence of a gravitational perturber and to acquire a con-verged value of the disk torque in the absence of the accretionheating, corresponding to the cold-protoplanet limit. For the κ BL -disk, this part of the simulation is similar to the experiments inLega et al. (2014). Subsequently, we continue the simulation foranother 30 P orb during which we let the planetary mass growwhile releasing the accretion heat into the gas disk accordingto Eqs. (7) and (8).Fig. 2 shows the temporal evolution of the torque exerted onthe protoplanet by the κ const -disk and κ BL -disk. During the phasewithout accretion heating, the torque converges to a stationaryvalue during 10 P orb . The torque value in the κ BL -disk is morepositive compared to the κ const -disk, which is because the steeperradial decline of the entropy in the κ BL -disk enhances the posi-tive entropy-driven part of the corotation torque (Paardekooper& Mellema 2006; Baruteau & Masset 2008).When the accretion heating is activated in the κ const -disk, apositive contribution is added to the torque, in accordance withthe results of Benítez-Llambay et al. (2015). The torque slightlyoscillates at first, but the amplitude of the oscillations decreasesin time and becomes negligible at t = P orb .On contrary, when the accretion heating is activated in the κ BL -disk, the outcome of the heating torque becomes very dif-ferent. Strong oscillations of the disk torque are excited almostimmediately and they do not vanish in time; instead, their ampli- Article number, page 4 of 19. Chrenko and M. Lambrechts: Oscillatory migration of accreting protoplanets driven by a 3D distortion of the gas flow cold non-luminous protoplanet ( κ const -disk simulation, t = P orb )perturbed density temperaturehot luminous protoplanet ( κ const -disk simulation, t = P orb )perturbed density temperature excess Fig. 3.
Hydrodynamic quantities in the midplane of the κ const -disk close to the protoplanet. Top row corresponds to the cold protoplanet rightbefore the accretion heating is initiated (at t = P orb ), bottom row shows the steady state of gas around the hot protoplanet (at t = P orb ).The figure is constructed as a Cartesian projection of the spherical grid. The density maps (right) display the perturbation ( ρ − ρ ) /ρ relative tothe equilibrium disk ( t = P orb ), the temperature map for the cold protoplanet (top right) shows the absolute values, and the temperature mapfor the hot protoplanet (bottom right) shows the excess with respect to the cold protoplanet (by subtracting T ( t = P orb ) from T ( t = P orb )).The position of the protoplanet is marked with the cross, the extent of its Hill sphere is bordered by the circle. The green curves (right) showthe topology of streamlines in the frame corotating with the protoplanet. In the inertial frame, the protoplanet would orbit counterclockwise. Thestreamlines outwards from its orbit thus depict the flow directed from y > y <
0, the inward streamlines are oriented in the opposite direction.A detailed view of the streamlines is provided in Fig. 5 where we also sort them according to their type. tude remains the same. The arithmetic mean of the torque mea-sured in the time interval 30–60 P orb is ¯ Γ (cid:39) − . × − a Ω ,implying that the torque is more negative compared to the situ-ation without accretion heating (the mean torque in the interval20–30 P orb is ¯ Γ (cid:39) − . × − a Ω ). The amplitude of the vari-ations with respect to the mean value is (cid:39)± . × − a Ω andthe oscillation period is (cid:39) . P orb .The torque oscillations are unexpected and we thus focusthroughout the rest of the paper on finding the physical mech-anism that excites them. The occurrence of oscillations suggeststhat the gas distribution around the protoplanet is changing dur-ing the simulation, as we can demonstrate using the integral ex-pression for the disk torque Γ = (cid:90) disk (cid:16) r p × F g (cid:17) ⊥ d V , (11)where r p is the protoplanet’s radius vector, F g is the gravitationalforce of a disk element, the vertical component of the cross prod-uct is considered and we integrate over the disk volume. Only anon-zero azimuthal component of F g can lead to a non-vanishing cross product in the integral, thus any oscillations of Γ must berelated to a variation of F g ,θ . In other words, there must be anazimuthal asymmetry in the gas distribution with respect to theprotoplanet for the torque to be non-zero and only a temporalredistribution of the asymmetry can cause a torque oscillation.Our strategy throughout the rest of Sect. 3 is the following:First, we focus on the κ const -disk simulation in Sect. 3.3. Al-though similar simulations were analysed by Benítez-Llambayet al. (2015), our aim is to focus on the 3D gas flow that hasnot been described yet. Our findings are then expanded for the κ BL -disk simulation in Sect. 3.4 where we relate the gas redis-tribution to the oscillatory behaviour of the torque. Sect. 3.5 isdevoted to identifying the key physical mechanisms that a ff ectthe gas flow. In Sect. 3.6, we vary the disk opacity gradient andwe study its impact on the torque oscillation. Finally, we relaxthe assumption of a fixed orbit and explore how the protoplanetmigrates in Sect. 3.7. In this section, the κ const -disk simulation is analysed. Article number, page 5 of 19 & A proofs: manuscript no. paper
Fig. 3 compares the midplane density and temperature distribu-tion around the cold and hot protoplanet. We display the state ofthe simulation at t =
30 (i.e. at the final stage of the phase with-out accretion heating) and 60 P orb (i.e. at the final stage of thephase with accretion heating). The latter represents the steadystate of gas around the accreting protoplanet and we verified thatsuch a distribution is achieved early (at t (cid:39) P orb ) and does notgreatly evolve since then.For the non-luminous protoplanet, the gas state is in agree-ment with Lega et al. (2014) (see their Fig. 10 for a compari-son). The density distribution is not spherically symmetric withrespect to the protoplanet but there are two patches of slightlyoverdense gas along the outflow from the Hill sphere known asthe cold fingers (explained in Sect. 1). The temperature drop in-side the fingers is clearly apparent from the top right panel ofFig. 3.When the protoplanet becomes luminous, the gas distri-bution is modified and the heating torque arises. FollowingBenítez-Llambay et al. (2015) and Masset (2017), one can imag-ine the response of the gas to the heating from the protoplanet asfollows: First, an underdense disturbance appears close to theprotoplanet, with a characteristic scale length given by the linearperturbation model of Masset (2017): λ c = (cid:114) χ q Ω p γ , (12)where χ is the thermal di ff usivity and q is a dimensionless mea-sure of the disk shear ( q = / ff usion interplay is indeed observed inFig. 3, where we find the typical two-lobed distribution of hotunderdense gas around the luminous protoplanet and the positiveboost of the total torque (Fig. 2) confirms that the outer lobe isslightly more pronounced. The bottom right panel reveals themagnitude and spatial extent of the temperature excess, as wellas its skewed shape in the direction of the disk shear.However, we make a new observation here concerning thestreamlines of the flow that are overlaid in the temperature maps.It is obvious that the hot perturbation significantly changes thetopology 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 passthe protoplanet and a new set of spiral-like retrograde stream-lines appears. It is known that vertical motions play an important role forthe 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; Popo-vas 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 thevertical plane intersecting the protoplanet’s location. When theprotoplanet is cold, there is a vertical stream of gas descendingtowards the protoplanet and escaping as an outflow through themidplane, in accordance with e.g. Lambrechts & Lega (2017).In the presence of accretion heating, however, the directionof the gas flow above the protoplanet is reverted. It forms anoutflowing column while the midplane flow becomes directedtowards the protoplanet. Two overturning cells appear on eachside of the vertical column (although it is important to point outthat no such cells are apparent in the full 3D flow which is dis-cussed later). We notice that the hot perturbation is not spheri-cally symmetric but rather elongated in the direction of the col-umn, indicating that the envelope is not hydrostatic.cold protoplanet ( κ const -disk simulation, t = P orb )hot protoplanet ( κ const -disk simulation, t = P orb ) Fig. 4.
Temperature map around the cold protoplanet at t = P orb (top) and temperature excess around the hot protoplanet at t = P orb (bottom) in the κ const -disk. Vertical plane (perpendicular to the disk mid-plane) is displayed. The green arrows show the vertical velocity vectorfield. So far, we described two new findings that were not incorpo-rated in the existing descriptions of the heating torque (Benítez-Llambay et al. 2015; Masset 2017): The distortion of the stream-line topology and the reversal of the vertical motions. The flowdirection is directly linked to the heating torque because it de-termines the redistribution of the hot gas by advection and thuscontributes to the shape of the underdense regions. Therefore wefocus on the streamline topology in this section.The streamlines are calculated using the explicit first orderEuler integrator and the trilinear interpolation of the velocityfield. The interpolation allows us to obtain the velocity vectorat an arbitrary location within the spherical grid. The size of theintegration step is chosen so that the length integrated during asingle propagation does not exceed 0 . v φ = Article number, page 6 of 19. Chrenko and M. Lambrechts: Oscillatory migration of accreting protoplanets driven by a 3D distortion of the gas flow cold protoplanet ( κ const -disk simulation, t = P orb ) -0.15-0.1-0.05 0 0.05 0.1 0.15 5.05 5.1 5.15 5.2 5.25 5.3 5.35 a z i m u t h θ [ r ad ] radius r [au] i nn e r ou t e r frontrear r ad i u s -0.05 0 0.05 a z i m u t h c o l a t i t ude hot protoplanet ( κ const -disk simulation, t = P orb ) -0.15-0.1-0.05 0 0.05 0.1 0.15 5.05 5.1 5.15 5.2 5.25 5.3 5.35 a z i m u t h θ [ r ad ] radius r [au] i nn e r ou t e r frontrear r ad i u s -0.05 0 0.05 a z i m u t h c o l a t i t ude Fig. 5.
Detailed streamline topology in the κ const -disk simulation. Top row corresponds to the cold protoplanet at t = P orb and the bottom row tothe hot protoplanet at t = P orb . Rectangular projection in the spherical coordinates is used to display the disk midplane near the protoplanet (left)and the actual 3D flow (right). The colour of the curves distinguish individual sets of streamlines: inner circulating (yellow), outer circulating (red),front horseshoe (light blue), rear horseshoe (purple) and other (green). The thick black lines highlight the critical circulating streamlines closestto the protoplanet. The black cross and the ellipse mark the protoplanet’s location and Hill sphere; the black arrows indicate the flow directionwith respect to the protoplanet. In 3D figures (right), the dark blue hemisphere corresponds to the Hill sphere above the midplane. Additionally,the insets in the corners of the 3D figures provide a close-up of the upstream outer circulating streamlines viewed from a slightly di ff erent angle.The endpoints indicate where the flow exits the depicted part of the space and highlight if the initially coplanar streamlines descend towards theprotoplanet (top) or rather rise vertically (bottom). We emphasise that the streamlines in the 3D figures are generated above the midplane and donot directly correspond to those in the 2D figures. a useful visualisation, we emphasise that they – by construction– carry zero information about the adjacent vertical flows. Forthe 3D projections, the streamlines are generated slightly abovethe midplane (at φ = π/ − .
005 rad) to take into account non-zero v φ .Fig. 5 shows 2D and 3D streamlines in the κ const -disk, againfor the exact same simulation stages as in Figs. 3 and 4. In theplots, we distinguish the following types of streamline: First,there are circulating streamlines that do not cross the proto-planet’s corotation. To imagine the direction of the relative mo-tion, we point out that the gas on inner circulating streamlines moves faster than the protoplanet whereas the gas on outer circu-lating streamlines lags behind the protoplanet. Second, there arehorseshoe streamlines that make a single U-turn and cross thecorotation once at each side of the protoplanet. Such streamlinesahead of the protoplanet’s orbital motion form the front horse-shoe region, those located behind the protoplanet belong to therear horseshoe region. To outline the separatrices between thehorseshoe and circulating regions, we highlight the critical innerand outer circulating streamlines that are located closest to theprotoplanet. Finally, some of our plots contain streamlines thatdo not fall in any of the aforementioned categories. Article number, page 7 of 19 & A proofs: manuscript no. paper hot protoplanet ( κ BL -disk simulation) Fig. 6.
Evolution of the perturbed midplane gas density in the κ BL -disk simulation. The corresponding simulation time t is given by labels.Individual snapshots represent the state of gas when the total torque acting on the protoplanet is minimal (left), maximal (third) and oscillating inbetween (second and right). The streamlines are overplotted for reference. The figure is also available as an online movie 1, showing the temporalevolution from t =
30 to 33 P orb . When the protoplanet is non-luminous, the 2D midplanestreamlines in Fig. 5 do not exhibit any unexpected features. Thestagnation point (X-point) of the flow is located within the Hillsphere, which is also intersected by both horseshoe and circu-lating streamlines. In 3D, we notice that upon making their U-turn, the horseshoe streamlines vertically descend towards themidplane, as already pointed out by Fung et al. (2015) or Lam-brechts & Lega (2017). The descend is also exhibited by someof the circulating streamlines closest to the protoplanet.For the hot protoplanet, we now obtain a clear picture of thestreamline distortion. In the midplane, the following changes ap-pear: – Circulating streamlines cross a smaller portion of the Hillsphere. When passing the protoplanet, they are bent towardsit (unlike near the non-luminous protoplanet where they arerather deflected away), which is especially apparent for thecritical circulating streamlines. – The classical horseshoe streamlines detach from the Hillsphere and make their U-turns at greater azimuthal separa-tions. – Part of the streamlines originating downstream the horseshoeregions is captured inside the Hill sphere where it rotatesaround 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 Hillsphere. – When the circulating streamlines pass the protoplanet andbecome perturbed, they are also uplifted. This behaviour isexactly the opposite to the cold-protoplanet situation.
We now return to the κ BL -disk simulation for which we discov-ered the strong oscillations of the heating torque (Fig. 2). Investigating the evolution of the gas density, we find that theposition and size of the underdense lobes never becomes station-ary, as shown in Fig. 6 (see also the online movie 1). The figureshows a sequence of snapshots corresponding to t = .
3, 31 . .
20 and 32 . P orb .The first panel depicts the state when the total torque reachesits first minimum during the beginning of the oscillatory phase.There is a dominant underdense lobe located ahead of the proto-planet while the rear lobe almost disappears. Such a distributioncan be easily related to the strong negative torque: The excessof the gas mass behind the protoplanet (and the paucity of massahead of it) leads to an azimuthal pull acting against the orbitalmotion, imposing a negative torque.When the torque is reversing from negative to positive (sec-ond panel), both lobes are similarly pronounced. The rear oneseems to be located closer to the protoplanet. The third panelcorresponds to the torque maximum. The rear lobe is dominantand thus the overabundance of the gas ahead of the protoplanetmakes the torque positive. The final panel shows the state whenthe oscillating torque is approximately half-way from positive tonegative. The gas distribution indeed looks like a counterpart tothe second panel since both lobes are again apparent but the frontone is now closer to the protoplanet.The gas redistribution is clearly related to the topology of thestreamlines and to the position of the spiral-like flow. We noticethat the centre of the captured streamlines undergoes retrograde(‘clockwise’) rotation around the protoplanet. One underdenselobe is always associated with this rotating flow. Apparently, theredistribution of the hot gas by advection tends to favour the lobewhich is intersected by the majority of the captured streamlinesat a given time. In the first panel of Fig. 6, for example, the hotgas is transported more e ffi ciently into the front lobe, creatinga strong front-rear asymmetry between the lobes. In the third Article number, page 8 of 19. Chrenko and M. Lambrechts: Oscillatory migration of accreting protoplanets driven by a 3D distortion of the gas flow hot protoplanet ( κ BL -disk simulation)(a) t = . P orb (b) t = . P orb (c) t = . P orb -0.15-0.1-0.05 0 0.05 0.1 0.15 5.05 5.1 5.15 5.2 5.25 5.3 5.35 a z i m u t h θ [ r ad ] radius r [au] i nn e r ou t e r frontrear -0.15-0.1-0.05 0 0.05 0.1 0.15 5.05 5.1 5.15 5.2 5.25 5.3 5.35 a z i m u t h θ [ r ad ] radius r [au] i nn e r ou t e r frontrear -0.15-0.1-0.05 0 0.05 0.1 0.15 5.05 5.1 5.15 5.2 5.25 5.3 5.35 a z i m u t h θ [ r ad ] radius r [au] i nn e r ou t e r frontrear (d) t = . P orb (e) t = P orb (f) t = . P orb -0.15-0.1-0.05 0 0.05 0.1 0.15 5.05 5.1 5.15 5.2 5.25 5.3 5.35 a z i m u t h θ [ r ad ] radius r [au] i nn e r ou t e r frontrear -0.15-0.1-0.05 0 0.05 0.1 0.15 5.05 5.1 5.15 5.2 5.25 5.3 5.35 a z i m u t h θ [ r ad ] radius r [au] i nn e r ou t e r frontrear -0.15-0.1-0.05 0 0.05 0.1 0.15 5.05 5.1 5.15 5.2 5.25 5.3 5.35 a z i m u t h θ [ r ad ] radius r [au] i nn e r ou t e r frontrear Fig. 7.
Midplane streamline topology in the κ BL -disk simulation. The panels are labelled by the simulation time. The individual types of streamlinesare the same as in Fig. 5. The sequence (a) to (f) represents the transition between the states corresponding to the minimum and maximum of thetorque, respectively (see Fig. 2 to relate the panels to the torque evolution). panel, the situation is exactly the opposite and the rear lobe ismore pronounced.We point out that the continuous variations of the hot lobesand their alternating dominance are unexpected features of theheating torque which was previously thought to be strictly posi-tive (Benítez-Llambay et al. 2015). We again explore the streamline topology and its changes relatedto the redistribution of the hot gas. Fig. 7 shows the midplanestreamlines near the protoplanet for a selection of simulationtimes between t = . . P orb . The former correspondsto the torque minimum, the latter to the torque maximum. Thetime intervals between the individual panels are not uniform butrather selected to highlight the most interesting transitions.The sequence reveals the following features: – In panel (a), the streamlines are similar to the steady-state κ const -disk case (bottom of Fig. 5) in several ways, mostlyin the detachment of the classical horseshoe streamlines andin the existence of the retrograde streamlines captured fromdownstream horseshoe regions. But there are also importantdi ff erences: There are none inner circulating streamlines in-tersecting the Hill sphere and moreover, a larger part of thecaptured streamlines originate in the rear horseshoe regionwhile the front horseshoe region is almost disconnected. – In panel (b), the front horseshoe region becomes entirely dis-connected. The captured streamlines originate exclusively inthe rear horseshoe region. – In panel (c), also the outer circulating streamlines stop cross-ing the Hill sphere. Some of the streamlines that were cap-tured in (b) now overshoot the protoplanet and make a U-turnahead of it. The centre around which the captured stream-lines enclose becomes shifted behind the protoplanet. – In panel (d), the front X-point moves closer to the Hill sphereand so do the front horseshoe streamlines. The rear horse-shoe region becomes radially narrower and the number ofU-turn streamlines overshooting the protoplanet diminishes. – In panel (e), the front horseshoe region reconnects with thecaptured streamlines. – In panel (f), the captured streamlines originate mostly inthe front horseshoe region while the rear horseshoe regionis evolving towards its disconnection, similarly to what wesaw for the front horseshoe region in panels (a) and (b). Thecentre of the captured streamlines moves inwards from theprotoplanet (and will continue to propagate ahead).Fig. 8 shows three selected snapshots of 3D streamlines( t = .
3, 31 .
75 and 32 . P orb ) when the oscillating torque is atits minimum (top), grows halfway towards the maximum (mid-dle) and reaches it (bottom). Clearly, the reshaping of the stream-line topology that we described for the midplane propagates in acomplicated way into the vertical direction as well: – In the first panel, the spiraling streamlines of the vertical col-umn are more or less centred above the protoplanet and as
Article number, page 9 of 19 & A proofs: manuscript no. paper they rise above the midplane they penetrate the majority ofthe Hill sphere. – In the second panel, we see the overshooting rear horseshoestreamlines making their U-turns within the Hill sphere. Thevertical column of captured streamlines is displaced to therear of the Hill sphere. As the captured streamlines spiral up,the column tilts towards the Hill sphere. The outer circulatingstreamlines are strongly uplifted towards colatitudes abovethe Hill sphere. – In the third panel, some uplifted outer circulating streamlinespenetrate into the vertical column and by this reconfigura-tion, the column reconnects with the front horseshoe regionwhile the rear one starts to disconnect.
In previous sections, we revealed that the gas heating from an ac-creting protoplanet changes the topology of the flow. Perturbedstreamlines have a tendency to bent towards the protoplanet andalso to rise vertically. If the perturbations become strong enough,the streamlines can form a vertical spiral. In this section, we in-vestigate the physical processes responsible for such a streamlinedistortion.In Sect. 3.5.1, we theorise that the streamline distortion isa result of vorticity perturbations arising because the vigorousaccretion heating renders the circumplanetary gas baroclinic. Aconfirmation is provided for the steady state of the κ const -disksimulation in Sects. 3.5.2 and 3.5.3. Finally, Sect. 3.5.4 demon-strates that the vertical temperature gradient above the accretingprotoplanet is superadiabatic and we also highlight di ff erencesbetween the κ const -disk and κ BL -disk. The spiral-like structure of the captured streamlines and thebending of nearby circulating streamlines suggests that the vor-ticity of the flow is modified when the protoplanet becomes hot.The vorticity can be defined via the relation ω a ≡ ω r + Ω ≡ ∇ × v + Ω , (13)where ω a is the absolute vorticity and ω r is the relative vorticityin the reference frame corotating with the protoplanet.Evolution of ω r is described by the vorticity equation in thecorotating frame (see Appendix C)D ω r D t = ( ω a · ∇ ) v − ω a ( ∇ · v ) + ∇ ρ × ∇ P ρ , (14)where D / D t denotes the Lagrangian derivative. In writing theequation, we neglected the e ff ects of viscous di ff usivity (largeReynolds number limit) but otherwise the equation is general.Regarding the right-hand side terms, the first one describesthe tendency of vortex tubes to become twisted due to veloc-ity field gradients. The second one characterises the stretchingor contraction of vortex tubes due to flow expansion or com-pression. These first two terms, usually called the twisting andstretching term, are only important if there is non-zero absolutevorticity already existing in the flow and they can cause its re-distribution.The remaining term on the right-hand side is the baroclinicterm. It vanishes in barotropic flows where the pressure anddensity gradients are always parallel, but since our model is hot protoplanet ( κ BL -disk simulation) t = . P orb r ad i u s -0.05 0 0.05 a z i m u t h c o l a t i t ude t = . P orb r ad i u s -0.05 0 0.05 a z i m u t h c o l a t i t ude t = . P orb r ad i u s -0.05 0 0.05 a z i m u t h c o l a t i t ude Fig. 8.
3D streamlines in the κ BL -disk simulation. Each snapshot is la-belled by the simulation time. The panels correspond to the minimum(top) and maximum (bottom) torque and to the state in between (mid-dle).Article number, page 10 of 19. Chrenko and M. Lambrechts: Oscillatory migration of accreting protoplanets driven by a 3D distortion of the gas flow cold non-luminous protoplanet ( κ const -disk simulation, t = P orb )azimuthal components radial components polar components -3.0-2.5-2.0-1.5-1.0-0.50.00.5 ω r , θ [ - c . u .] -3.0-2.5-2.0-1.5-1.0-0.50.00.51.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 s ou r c e s [ - c . u .] time [c.u.]baroclinicstretchingtwisting -0.4-0.20.00.20.40.60.81.01.21.41.61.8 ω r , r [ - c . u .] -1.0-0.50.00.51.01.52.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 s ou r c e s [ - c . u .] time [c.u.] 10.511.011.512.012.513.013.514.0 ω r , φ [ - c . u .] -1.5-1.0-0.50.00.51.01.52.02.53.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 s ou r c e s [ - c . u .] time [c.u.] hot luminous protoplanet ( κ const -disk simulation, t = P orb )azimuthal components radial components polar components -1.0-0.50.00.51.01.52.02.53.03.5 ω r , θ [ - c . u .] -1.0-0.50.00.51.01.52.02.53.03.5 60.0 80.0 100.0 120.0 140.0 160.0 s ou r c e s [ - c . u .] time [c.u.]baroclinicstretchingtwisting -0.3-0.2-0.10.00.10.20.30.4 ω r , r [ - c . u .] -1.5-1.0-0.50.00.51.01.5 60.0 80.0 100.0 120.0 140.0 160.0 s ou r c e s [ - c . u .] time [c.u.] 10.010.511.011.512.012.513.013.514.0 ω r , φ [ - c . u .] -5.0-4.0-3.0-2.0-1.00.01.02.0 60.0 80.0 100.0 120.0 140.0 160.0 s ou r c e s [ - c . u .] time [c.u.] Fig. 9.
Evolution of the relative vorticity (first and third row) and balance of the vorticity source terms (second and fourth row) along a singlestreamline in the κ const -disk simulation. The streamlines for these measurements are chosen from the 3D sets displayed in Fig. 5. For the hotprotoplanet (bottom two rows), the extremal outer circulating streamline is selected and for the cold protoplanet (top two rows), we select anouter circulating streamline with a comparable Hill sphere crossing time. Azimuthal (left), radial (middle) and polar (right) components of thevorticity (black solid curve), baroclinic term (blue solid curve), stretching term (orange dashed curve) and twisting term (magenta dotted curve)are displayed using scaled code units. The grey rectangle marks the Hill sphere crossing. We point out that the source terms represent the rate ofchange of the vorticity and also that the vertical range is not kept fixed among individual panels. not barotropic, ∇ ρ and ∇ P can be misaligned, leading to vor-ticity production or destruction since their cross product canbe non-zero. Because ρ near the accreting protoplanet exhibitsasymmetric perturbations while P remains roughly sphericallysymmetric, one can expect non-zero baroclinic perturbations toarise. The respective non-zero vorticity then enhances circula-tion around a given point of the continuum, twisting the stream-lines with respect to the situation unperturbed by accretion heat-ing. It is not a priori evident which source term is the most importantfor the vorticity evolution in our simulations. In Fig. 9, we studythe variations of ω r and the source terms of Eq. (14) along a sin-gle outer circulating streamline. We compare the situation nearthe cold and hot protoplanet in the κ const -disk.Downstream, before the streamline encounters the proto-planet, the situation is similar for the compared cases: the az-imuthal and radial vorticity components ω r ,θ and ω r , r are negli-gible, while the polar component ω r ,φ is non-zero and positive.The positive value of ω r ,φ corresponds to the inherent vorticityin a flow with the Keplerian shear: Taking v r = v φ = v θ = √ GM / r − Ω r , one obtains ω r ,φ = − . Ω K ( r ) + Ω . At thesame time, the source terms are zero because far from the pro-toplanet there are no strong velocity gradients, no compressionand ∇ ρ and ∇ P are aligned.To relate the vorticity variation with the streamline distortionclose to the protoplanet, let us make a thought experiment, utilis-ing 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. Imaginean observer moving along the critical outer circulating stream-line, corresponding to the outer thick black curve in the bottomright panel of Fig. 5. The observer moves with the flow, pre-dominantly in the − θ direction. For this experiment, we dub thedirections θ , − θ , r , − r , φ and − φ as behind, ahead, outwards, in-wards, down and up, respectively. Considering only the stream-lines of Fig. 5 originating at r > r p , they initially define a planeat constant φ = π/ − .
005 rad and the observer propagatesthrough that plane. When the flow reaches the Hill sphere, theobserver studies the instantaneous displacement of nearby gasparcels which will correspond to the deformation of the surfacedefined by neighbouring streamlines.According to bottom panels of Fig. 9, the observer measures ω r ,θ > ω r ,θ points againstthe direction of the observer’s motion, forcing the circulationin the local ( r , φ ) plane. The nearby gas parcels must obey thefollowing right-hand rule: When the thumb points in the direc-tion of ω r ,θ , wrapping fingers determine the direction of circula-tion. From the observer’s point of view, an outer gas parcel willbe falling downwards to the midplane and an inner gas parcelwill be rising upwards. This is exactly in accordance with Fig. 5where streamlines passing the protoplanet are uplifted from the φ = π/ − .
005 rad plane and the kick gets stronger with de-creasing separation to the protoplanet. ω r , r only slightly oscillates during the Hill sphere crossingbut becomes positive (although relatively small) upstream. Us-ing again the same considerations as above, ω r , r points outwardsfrom the observer after crossing the Hill sphere, promoting cir- Article number, page 11 of 19 & A proofs: manuscript no. paper culation in the ( θ, φ ) plane. Using the right-hand rule, a gas par-cel ahead of the observer will be falling downwards and a gasparcel behind the observer will be rising upwards. In Fig. 5, thisis reflected by the streamline topology when the red streamlinesrising after the Hill sphere passage suddenly start to fall backtowards the midplane.As for the remaining vorticity component ω r ,φ , it remainspositive during the Hill sphere crossing, but it acquires a positiveboost at first and a more prominent negative perturbation after-wards. The later diminishes the circulation related to shear in the( r , θ ) plane. This is only possible if gas parcels near the observerbecome displaced towards trajectories with smaller shear veloci-ties. In Fig. 5, the streamline topology indeed exhibits such a be-haviour because when the red streamlines (and similarly greenand yellow ones) pass the protoplanet, they are being bent to-wards it.Looking at the source terms, it is obvious that the baroclinicterm is responsible for perturbing ω r ,θ and ω r ,φ , while counteract-ing the twisting term contributing to ω r , r . The importance of thebaroclinic term for the flow approaching the hot protoplanet isthus confirmed. Moreover, the perturbation is indeed 3D as eachof the studied components is important for the resulting stream-line topology.Comparing the hot-protoplanet case to the cold-protoplanetcase, we notice that the evolution of ω r ,θ and ω r ,φ is roughly anti-symmetric, as well as the evolution of the baroclinic source term.This is consistent with our finding that streamlines near the coldprotoplanet are distorted in the opposite manner (they fall to-wards the midplane and deflect away from the protoplanet). Inother words, the baroclinic behaviour of the protoplanet’s vicin-ity is reverted between the cold- and hot-protoplanet case. Although we do not repeat the vorticity analysis for the remain-ing sets of streamlines, it is clear that features of the streamlinedistortion can be explained by the baroclinic generation of thevorticity. To further support this claim, we compare in Fig. 10the map of the baroclinic term near the cold and hot protoplanetin the κ const -disk. We plot the polar component of the baroclinicterm ( ∇ ρ × ∇ P ) φ /ρ in the midplane. Since the midplane flowis e ff ectively 2D (because v φ = ω r is non-zero and thus also the polar componentof the baroclinic term is the most important one.Fig. 10 reveals that the gas is baroclinic near both the coldand hot protoplanet, but for the latter case, the map becomesapproximately antisymmetric compared to the cold-protoplanetcase, slightly rotated in the retrograde sense and the baroclinicregion is more extended.The existence of the baroclinic region can be explained usingthe isocontours of constant volume density and isobars. For abarotropic gas ( ∇ ρ (cid:107) ∇ P ), any nearby isocontours of constant ρ and P should have the same shape. Wherever the isocontoursdepart from one another, it means that the local gradients of ρ and P are misaligned and that the baroclinic term is non-zero.Looking at the cold-protoplanet case, the isobars and con-tours of constant density are nearly spherically symmetric. How-ever, we notice that each density isocontour exhibits two bumpsand appears to be stretched in the direction where the gas out-flows from the Hill sphere. Clearly, these bumps are associatedwith the cold-finger perturbation that appears in the same loca-tion (see Fig. 3). Since the cold fingers are filled with overdensegas, they perturb the local density gradient, making it to point to- wards them. At the same time, the isobars remain approximatelyspherically symmetric.When the protoplanet is hot, the cold fingers are replacedwith hot underdense perturbations. Therefore, the local densitygradient tends to point away from them. We can see that thisis indeed true because the density isocontours do not exhibitbumps, but rather concavities across the overheated region (com-pare with Fig. 3). Clearly, this is the reason why the baroclinicmap for the hot protoplanet looks reverted compared to the coldprotoplanet.Summarising these findings, we saw that the thermal pertur-bations associated either with the cold-finger e ff ect or the heatingtorque make the circumplanetary region baroclinic. But the influ-ence on the vorticity evolution is the opposite when comparingthe cold- and hot-protoplanet case.cold protoplanet ( κ const -disk simulation, t = P orb )hot protoplanet ( κ const -disk simulation, t = P orb ) Fig. 10.
Maps of the φ -component of the baroclinic term in the κ const -disk simulation. The purple isocontours depict several levels of the con-stant volume density and the green isocontours correspond to the iso-bars. The levels of the contours are kept fixed between the panels. Although the baroclinic distortion of the gas flow is a robustmechanism, it does not provide a simple explanation why the κ BL -disk simulation exhibits gas instability whereas the κ const -disk simulation remains stable. We now explore the vertical sta-bility of both disks against vertical convection, considering onlythe hot-protoplanet limit. Article number, page 12 of 19. Chrenko and M. Lambrechts: Oscillatory migration of accreting protoplanets driven by a 3D distortion of the gas flow
To do so, we employ the Schwarzschild criterion |∇ rad ,φ | > |∇ ad ,φ | ⇔ |∇ rad ,φ ||∇ ad ,φ | − > , (15)where the subscript ‘rad’ denotes the vertical temperature gra-dient found in our simulations and the subscript ‘ad’ denotes thetemperature gradient that an adiabatic gas would establish. The ∇ symbol stands for the logarithmic gradient d log T / d log P , yield-ing |∇ ad ,φ | = ( γ − /γ for the adiabatic case.The Schwarzschild criterion is not necessarily a universalway to determine if the disk is unstable to convection. The rea-son is that convective destabilisations in protoplanetary disks areopposed by di ff usive e ff ects (e.g. Held & Latter 2018) and shearmotions (Rüdiger et al. 2002). However, to our knowledge thereare no convective criteria that would take into account the diskperturbation by the protoplanet, we thus use the Schwarzschildcriterion for its simplicity, keeping the limitations in mind.Fig. 11 shows the vertical velocity field and balance of theSchwarzschild criterion in the vertical direction of the κ const -diskand κ BL -disk. Each panel shows a di ff erent simulation time andvertical plane. For the κ const -disk, we display t = P orb andthe vertical plane intersecting the protoplanet’s location, whereasfor the κ BL -disk, each plane approximately intersects the centrearound which the captured streamlines of Fig. 6 circulate at thegiven simulation time ( t = .
3, 31 .
75, 32 . . P orb ). Inother words, we choose the vertical planes where we expect themost prominent vertical outflow.The first thing we point out is that the background di ff ersbetween the studied disks. The background of the κ BL -disk isslightly superadiabatic, contrary to the κ const -disk. The di ff erencearises as a result of the opacity laws. As derived by Lin & Pa-paloizou (1980) and Ruden & Pollack (1991), one can make aqualitative estimate for optically thick regions unperturbed bythe protoplanet (cid:32) |∇ rad ,φ ||∇ ad ,φ | − (cid:33) background = / (4 − β )( γ − / ( γ ) − , (16)where β is the power-law index of the κ ∝ T β dependence. Inthe κ const -disk, β = |∇ rad ,φ | / |∇ ad ,φ | − (cid:39) − .
17. In the κ BL -disk, the Bell & Lin (1994) opacity law in the given temper-ature range corresponds to water-ice grains and exhibits β = |∇ rad ,φ | / |∇ ad ,φ | − (cid:39) .
67. These estimated values arein a good agreement with the background values of Fig. 11.The background itself is not convective because vertical con-vection is usually not self-sustainable unless there is a strongheat deposition within the disk (e.g. Cabot 1996; Stone & Balbus1996; Klahr et al. 1999; Lesur & Ogilvie 2010) which, however,can be provided by the hot accreting protoplanet in our case. In-deed, looking at the κ const -disk in Fig. 11 (top), the region wherewe previously identified the most significant temperature excessdue to planetary luminosity (compare with Fig. 4) is superadia-batic and the corresponding vertical outflow can be consideredconvective.In the κ BL -disk simulation, the temperature gradient departsfrom the adiabatic one even more and additionally, the excess isno longer centred above the protoplanet itself but rather spans itsvicinity. This can be seen from the varying azimuthal coordinateof the displayed vertical planes and also from the radial o ff set ofthe highly superadiabatic region in panels 3 and 5. We use the colatitude φ to study the vertical gradients because thecurvature of spherical coordinates near the midplane does not signifi-cantly depart from the true vertical direction. κ const -disk, hot protoplanet ( t = P orb ) κ BL -disk, hot protoplanet ( t = .
3, 31 .
75, 32 . . P orb ) Fig. 11.
Balance of the Schwarzschild criterion in the vertical planes ofthe κ const -disk (top) and the κ BL -disk (remaining panels). The individ-ual vertical planes are chosen to track the most prominent vertical flowat the given simulation time t (see the main text) and their azimuthalseparations from the protoplanet location are given by panel labels (us-ing multiples of the azimuthal span ∆ θ H of the Hill sphere radius). Thecolour maps evaluate Eq. (15). Positive values indicate superadiabaticvertical temperature gradients. The vertical velocity vector field is over-laid in the plots. The half-circles mark the overlap of a given plane withthe Hill sphere (the planes in panels 3 and 5 do not overlap with the Hillsphere). Article number, page 13 of 19 & A proofs: manuscript no. paper
We summarise the section by speculating that the κ BL -disksimulation becomes destabilised because the hot disturbance cre-ated by the accreting protoplanet is subject to vertical buoyantforces acting over a more extended region (compared to the κ const -disk). The reason is that the hot disturbance is imposedover an already-superadiabatic background. The uplift of thematerial cannot be compensated for in a stationary manner andeventually, the vertical outflow becomes o ff set with respect tothe protoplanet and starts to change its position in a cyclic man-ner. However, such a description is rather qualitative and preciseconditions for triggering the instability should be explored in fu-ture works. We now perform a partial exploration of the parametric spaceby varying the opacity gradient within the disk. The aim istwofold: First, we would like to support the claim of the previousSect. 3.5.4 about the importance of the vertical stratification forthe torque oscillations. Second, it is desirable to show that the ap-pearance of torque oscillations in the κ BL -disk is not coincidentaland that it can be recovered for a wider range of parameters.We construct six additional disk models with artificial opac-ity laws that i) conserve the opacity value at the protoplanet lo-cation (1 .
11 cm g − ) and ii) lead to opacity gradients which areintermediate between the κ const - and κ BL -disks. The latter prop-erty accounts for the highly unconstrained size distribution ofsolid particles in protoplanetary disks which manifests itself in alarge parametric freedom of the power-law slope of the opacityprofile (e.g. Piso et al. 2015).The first additional set of disks utilises the opacity law whichwe dub T -dependent: κ (cid:16) ¯ T ( r , φ ) (cid:17) = κ ¯ T β . (17)Similarly to the Bell & Lin (1994) opacity in the water-iceregime, it is exclusively a function of temperature. The tempera-ture ¯ T is again azimuthally averaged to disentangle the influenceof global opacity gradients (which we focus on) from those re-lated to accretion heating of the protoplanet. We examine thevalues β = .
5, 1 and 0 . β = κ const -disk) and 2 ( κ BL -disk). The constant of proportionality κ is always chosen to recover κ = .
11 cm g − at r = a p and φ = π/ κ (cid:39) . .
014 and0 .
122 cm g − for the respective values of β .The opacity law used for the second additional set of disks,referred to as r -dependent, is κ ( r ) = . (cid:32) ra p (cid:33) − δ cm g − , (18)and we choose δ =
3, 2 and 1 (because δ = κ BL -disk). The motivation for choosing thispurely radially dependent opacity law is to distinguish betweene ff ects caused by radial and vertical opacity gradients. The latterwill not appear when κ = κ ( r ).The opacity profiles of these disks in radiative equilibriumare summarised in Fig. 12 which reveals that all radial opac-ity gradients (top two panels) are indeed shallower compared tothe κ BL -disk. However, all disks with r -dependent opacities havezero vertical opacity gradient by construction (bottom panel),unlike disks with T -dependent opacities which vertically de-crease.The torque measurements are done the same way as in ourprevious experiments and the results are given in Fig. 13. For disks with T -dependent opacities (top), we find that the torqueoscillations appear for all investigated values of β . The oscilla-tion amplitude, however, linearly decreases with β and the periodbecomes slightly shorter as well. In case of r -dependent opaci-ties (bottom), the torque evolution does not strongly depend on δ and exhibits marginal and vanishing oscillations, as in the κ const -disk case.Since the only qualitative di ff erence between the disks with T -dependent and r -dependent opacities is in the vertical opacitygradient, our torque measurements confirm the importance of thevertical structure for the torque oscillations. The torque oscillateswherever the vertical opacity profile favours superadiabatic ver-tical stratification and the oscillation amplitude scales with thestrength of the opacity-temperature coupling. In the absence ofthe vertical opacity gradient, the oscillations are not established. opa c i t y κ [ c m g - ] protoplanetlocation T -dependentopacities κ BL κ ~ T κ ~ T κ ~ T r [au] r -dependentopacities κ BL κ ~ r -3 κ ~ r -2 κ ~ r -1 opa c i t y κ [ c m g - ] colatitude φ [rad/ π ] κ ~ r - δ κ ~ T κ ~ T κ ~ T Vertical profiles atthe protoplanet location
Fig. 12.
Radial (top two panels) and vertical (bottom) opacity profiles inequilibrium disks which we use to study the torque dependence on theopacity gradients. In top panels, the individual cases are distinguishedby colour and labelled in the legend. The profile of the κ BL -disk (solidgrey curve) is plotted for comparison. The bottom panel correspondsto the protoplanet location and demonstrates that only the disks with T -dependent opacities develop a vertical opacity gradient (which is notallowed for r -dependent opacities by construction). The previous simulations were conducted with the assumptionof the fixed orbit of the protoplanet and the static torque wasexamined. Here we explore whether our findings can be readilyapplied to a dynamical case when the protoplanet is allowed toradially migrate and its semimajor axis evolves. In this section,we focus only on the κ BL -disk in which we found the flow insta-bility.Starting from t = P orb , we release the protoplanet andrun the simulation until t = P orb . Fig. 14 compares the ob-tained dynamical torque with the previous result of our static Article number, page 14 of 19. Chrenko and M. Lambrechts: Oscillatory migration of accreting protoplanets driven by a 3D distortion of the gas flow -4e-05-3e-05-2e-05-1e-05 0 1e-05 2e-05 3e-05 t o t a l t o r que Γ [ a J Ω J ] Hot protoplanetCold protoplanet T -dependentopacities κ BL κ ~ T κ ~ T κ ~ T -4e-05-3e-05-2e-05-1e-05 0 1e-05 2e-05 3e-0520 25 30 35 40 45 50 t o t a l t o r que Γ [ a J Ω J ] no. of orbits [ P orb ] r -dependentopacities κ BL κ ~ r -3 κ ~ r -2 κ ~ r -1 Fig. 13.
Torque evolution in disks described in Fig. 12. The top panelcorresponds to disks with T -dependent opacities and the bottom panelto those with r -dependent opacities. The individual cases are distin-guished by colour and labelled in the legend. The evolution from the κ BL -disk simulation (solid grey curve) is given for reference. In the toppanel, the torque amplitude diminishes with the power-law index of theopacity law, yet the oscillations appear in all studied cases. In the bot-tom panel, we find that oscillations are rapidly damped. experiments. It is obvious that the oscillating character of thetorque is retained and therefore the instability of the flow op-erates near a moving protoplanet as well. There are di ff erencesboth in the amplitude and phase of the torque oscillations, butthe mean value of the torque over the simulated period of timeis ¯ Γ (cid:39) − . × − a Ω . Although this value is slightly morepositive than the static heating torque ( ¯ Γ (cid:39) − . × − a Ω ),it is almost the same as the torque acting on the cold protoplanet( ¯ Γ (cid:39) − . × − a Ω ). We thus confirm that in the κ BL -disk, theheating torque does not add any considerable positive contribu-tion to the mean torque but makes it strongly oscillating instead.Bottom panel of Fig. 14 shows the actual evolution of theprotoplanet’s semimajor axis. On average, the protoplanet slowlymigrates inwards but this drift is not smooth. The protoplanetexhibits fast periodic inward and outward excursions on an or-bital time scale. The migration rate of these individual excur-sions (not to be confused with the mean migration rate statedabove) is ˙ a ∼ (10 − au) / P orb .Finally, we note that the mean migration rate is not constant,which also corresponds to the varying o ff set of the dynamicaltorque with respect to the static torque in Fig. 14. Probably, theunstable gas distribution around the protoplanet is further af-fected by the protoplanet’s radial drift.
4. Discussion
This section discusses the applicability of our results, drawslinks to some previous studies and speculates about possible im- -6e-05-4e-05-2e-05 0 2e-05 4e-05 30 35 40 45 50 55 60 t o t a l t o r que Γ [ a J Ω J ] no. of orbits [ P orb ]Hot protoplanetdynamicalstatic s e m i m a j o r a x i s a [ au ] no. of orbits [ P orb ] κ BL , dynamical case Fig. 14.
Top: Comparison of the static (solid black curve) and dynami-cal torque (dashed blue curve) acting on the hot protoplanet in the κ BL -disk. Bottom: Evolution of the semimajor axis in the κ BL -disk when theprotoplanet is allowed to migrate. The migration is inward and oscilla-tory. plications for planet formation. Additionally, we outline possi-bilities for future work. Although most of them are beyond thescope of this paper, we at least present several additional sim-ulations 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 evolutionof torque oscillations). Our results revealed unexpected perturbations of the gas flow inthe vicinity of the protoplanet undergoing strong accretion heat-ing, one of them being the vertical outflow. The resolution of oursimulations was originally tailored for studying the torque andmay lack some accuracy close to the protoplanet where the out-flow occurs. Follow-up studies should thus utilise simulationswith high resolution of the protoplanetary envelope, similarlyto 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 κ andluminosity L for a given planetary mass M p exists for whichthe vertical outflow is triggered, similarly to the dependenciesof the heating torque explored by Benítez-Llambay et al. (2015).Popovas et al. (2018, 2019) studied the stability of the circum-planetary envelope during pebble accretion and found that thegas within the Bondi sphere exhibits 3D convective motions, as-suming L (cid:39) . × erg s − , M p = . M ⊕ and κ = g − .Lambrechts & Lega (2017) also explored a set of ( L , κ , M p ) pa-rameters and its impact on the structure of the circumplanetary Greater number of parameters were discussed in Popovas et al.(2018, 2019) but we quote those closest to this paper.Article number, page 15 of 19 & A proofs: manuscript no. paper envelope. They found the inner region of the envelope to de-part from hydrostatic equilibrium when the luminosity exceeds L = erg s − around a M p = M ⊕ core within a κ = g − environment, but they did not identify any vertical outflow. Theoutflow in our simulations appears for higher L and smaller M p compared to Lambrechts & Lega (2017), we can therefore as-sume that our parameters cross the critical ones.Regarding the baroclinic perturbations, although they areknown 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), theyhave been rarely considered in relation to hot protoplanets. Forexample, Owen & Kollmeier (2017) claim that that hot proto-planets can excite large-scale baroclinic vortices but we do notidentify any of those in our simulations. Instead, we find baro-clinic perturbations to be responsible for 3D distortion of the gasflow near the protoplanet. Although the heating torque has been previously thought to bestrictly positive and also e ffi cient in high-opacity locations ofprotoplanetary disks, our paper shows that it can exhibit morecomplicated behaviour if the temperature dependence of the diskopacity is taken into account. We identified an oscillatory modeof the heating torque in the disk region with κ ∝ T and wedemonstrated that it can operate even for shallower dependen-cies such as κ ∝ T . , albeit with a decreased amplitude. This be-haviour resembles the nature of baroclinic and convective diskinstabilities which usually operate in the most opaque regionswith the steepest entropy gradients but become less e ff ectiveelsewhere (e.g. Pfeil & Klahr 2019).It is worth noting that our simulations neglected the e ff ectof stellar irradiation. Stellar-irradiated disks tend to have verti-cal temperature profiles that are closer to being isothermal (e.g.Flock et al. 2013), unlike disk models used in this work whichhave rather adiabatic or slightly superadiabatic vertical tempera-ture gradients. On the other hand, even the irradiated disks oftencontain shadowed regions protected against stellar irradiation,for example behind the pu ff ed-up inner rim (e.g. Dullemondet al. 2001) or between the viscously heated inner disk and aflared irradiated outer disk (e.g. Bitsch et al. 2013). In such re-gions, the oscillatory migration could still operate.For the aforementioned reasons, it is likely that transitionzones might exist in protoplanetary disks, separating regionswhere protoplanets migrate under the influence of the standardpositive heating torque and where they undergo oscillatory mi-gration. Migration at the edges of such zones could be conver-gent, leading to a pileup of protoplanets.
5. Conclusions
By means of 3D radiation-hydrodynamic simulations, we inves-tigated the heating torque (Benítez-Llambay et al. 2015) actingon a luminous 3 M ⊕ protoplanet heated by accretion of solids.The aim was to compare the torque evolution and physics in adisk with non-uniform opacities (Bell & Lin 1994) with the out-come of a constant-opacity simulation.We discovered that the gas flow near the protoplanet is per-turbed by two mechanisms:1. The gas advected past the protoplanet becomes hot and un-derdense. Consequently, a misalignment is created betweenthe gradients of density and pressure within the protoplanet’s Hill sphere. The baroclinic term of the vorticity equation( ∼∇ ρ × ∇ P ) then becomes non-zero and modifies the vor-ticity of the flow.2. The e ffi cient heat deposition in the midplane makes the ver-tical temperature gradient superadiabatic, thus positively en-hancing vertical gas displacements.The streamline topology exhibits a complex 3D distortion. Themost important feature are spiral-like streamlines rising verti-cally above the hot protoplanet, forming an outflow column ofgas escaping the Hill sphere.In the constant-opacity disk, the vertical outflow is cen-tralised above the protoplanet, it temporarily captures stream-lines from both horseshoe regions and such a state is found tobe stationary over the simulation time scale. The distribution ofthe hot gas then remains in accordance with findings of Benítez-Llambay et al. (2015), having a two-lobed structure, and so doesthe resulting positive heating torque.In the disk with non-uniform opacity κ ∝ T (typically out-side the water-ice line), we find the superadiabatic temperaturegradient to be steeper and the distorted gas flow to be unstable.The vertical spiral flow becomes o ff set with respect to the pro-toplanet and periodically changes its position, spanning the edgeof the Hill sphere in a retrograde fashion. Its motion is followedby the underdense gas and the resulting heating torque stronglyoscillates in time. The interplay can be characterised by the fol-lowing sequence:1. A stage when most of the captured streamlines originate inthe rear horseshoe region and their spiral-like structure is o ff -set ahead of the protoplanet. Therefore the hot gas cumu-lates ahead of the protoplanet, a dominant underdense lobeis formed there and the torque becomes negative, reachingthe minimum of its oscillation.2. A stage when the front horseshoe region becomes com-pletely isolated from the captured streamlines. Some of therear horseshoe streamlines start to overshoot the protoplanetand make U-turns ahead of it. At the same time, the spiral-like structure recedes behind the protoplanet. The lobe fromstage 1 starts to decay while a rear lobe starts to grow and thetorque changes from negative to positive.3. Antisymmetric situation to stage 1, when most of the cap-tured streamlines originate in the front horseshoe region, thedominant lobe trails the protoplanet and the torque is posi-tive, reaching the maximum of its oscillation.4. Antisymmetric situation to stage 2, when the torque de-creases from positive to negative and the cycle repeats.Such an advective redistribution of the hot underdense gas is sus-tained over the simulation time scale.We also studied the dependence of the torque oscillationson the opacity gradient and found that they can appear even for κ ∝ T . , although their amplitude linearly decreases with thepower-law slope of the κ ( T ) dependence. We also demonstratedthat the oscillations would vanish in a disk with zero verticalopacity gradient.If the protoplanet is allowed to migrate, its mean migrationrate is nearly una ff ected but the radial drift is not smooth. It israther oscillatory, consisting of brief inward and outward radialexcursions with a characteristic rate ˙ a ∼ − au yr − .We discussed possible implications of the oscillating heat-ing torque for planet formation and pointed out that it can af-fect the global evolution of hot migrating low-mass protoplanets.During their migration through disk regions with varying opac-ities it might be possible for the protoplanets to switch between Article number, page 16 of 19. Chrenko and M. Lambrechts: Oscillatory migration of accreting protoplanets driven by a 3D distortion of the gas flow the standard positive heating torque of Benítez-Llambay et al.(2015) and the positive / negative torque oscillations discoveredin this paper. Acknowledgements.
The work of OC has been supported by Charles Univer-sity (research program no. UNCE / SCI / fargoca for comparison simulations. The authors are grateful toan anonymous referee whose insightful ideas helped to improve the manuscript. References
Barge, P., Richard, S., & Le Dizès, S. 2016, A&A, 592, A136Baruteau, C., Fromang, S., Nelson, R. P., & Masset, F. 2011, A&A, 533, A84Baruteau, C. & Masset, F. 2008, ApJ, 672, 1054Bell, K. R. & Lin, D. N. C. 1994, ApJ, 427, 987Benítez-Llambay, P., Masset, F., Koenigsberger, G., & Szulágyi, J. 2015, Nature,520, 63Benítez-Llambay, P. & Masset, F. S. 2016, ApJS, 223, 11Bitsch, B., Crida, A., Morbidelli, A., Kley, W., & Dobbs-Dixon, I. 2013, A&A,549, A124Brasser, R., Matsumura, S., Muto, T., & Ida, S. 2018, ApJ, 864, L8Brož, M., Chrenko, O., Nesvorný, D., & Lambrechts, M. 2018, A&A, 620, A157Cabot, W. 1996, ApJ, 465, 874Chrenko, O., Brož, M., & Lambrechts, M. 2017, A&A, 606, A114Chrenko, O., Brož, M., & Nesvorný, D. 2018, ApJ, 868, 145Cimerman, N. P., Kuiper, R., & Ormel, C. W. 2017, MNRAS, 471, 4662Coleman, G. A. L. & Nelson, R. P. 2016, MNRAS, 457, 2480Cresswell, P., Dirksen, G., Kley, W., & Nelson, R. P. 2007, A&A, 473, 329Crida, A., Morbidelli, A., & Masset, F. 2006, Icarus, 181, 587de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529Dullemond, C. P., Dominik, C., & Natta, A. 2001, ApJ, 560, 957Eklund, H. & Masset, F. S. 2017, MNRAS, 469, 206Flock, M., Fromang, S., González, M., & Commerçon, B. 2013, A&A, 560, A43Fung, J., Artymowicz, P., & Wu, Y. 2015, ApJ, 811, 101Fung, J. & Chiang, E. 2017, ApJ, 839, 100Goldreich, P. & Tremaine, S. 1979, ApJ, 233, 857Held, L. E. & Latter, H. N. 2018, MNRAS, 480, 4797Klahr, H. & Kley, W. 2006, A&A, 445, 747Klahr, H. H. & Bodenheimer, P. 2003, ApJ, 582, 869Klahr, H. H., Henning, T., & Kley, W. 1999, ApJ, 514, 325Kley, W. 1989, A&A, 208, 98Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971Kley, W. & Crida, A. 2008, A&A, 487, L9Korycansky, D. G. & Pollack, J. B. 1993, Icarus, 102, 150Kretke, K. A. & Lin, D. N. C. 2012, ApJ, 755, 74Kurokawa, H. & Tanigawa, T. 2018, MNRAS, 479, 635Lambrechts, M. & Johansen, A. 2012, A&A, 544, A32Lambrechts, M. & Johansen, A. 2014, A&A, 572, A107Lambrechts, M. & Lega, E. 2017, A&A, 606, A146Lega, E., Crida, A., Bitsch, B., & Morbidelli, A. 2014, MNRAS, 440, 683Lega, E., Morbidelli, A., Bitsch, B., Crida, A., & Szulágyi, J. 2015, MNRAS,452, 1717Lesur, G. & Ogilvie, G. I. 2010, MNRAS, 404, L64Lesur, G. & Papaloizou, J. C. B. 2010, A&A, 513, A60Levermore, C. D. & Pomraning, G. C. 1981, ApJ, 248, 321Li, H., Lubow, S. H., Li, S., & Lin, D. N. C. 2009, ApJ, 690, L52Lin, D. N. C. & Papaloizou, J. 1980, MNRAS, 191, 37Masset, F. S. 2002, A&A, 387, 605Masset, F. S. 2017, MNRAS, 472, 4204Masset, F. S. & Casoli, J. 2009, ApJ, 703, 857Masset, F. S. & Velasco Romero, D. A. 2017, MNRAS, 465, 3175McNally, C. P., Nelson, R. P., Paardekooper, S.-J., & Benítez-Llambay, P. 2019,MNRAS, 484, 728Mihalas, D. & Weibel Mihalas, B. 1984, Foundations of radiation hydrodynam-ics (Oxford University Press, New York)Miranda, R. & Lai, D. 2018, MNRAS, 473, 5267Morbidelli, A., Crida, A., Masset, F., & Nelson, R. P. 2008, A&A, 478, 929Ormel, C. W. & Klahr, H. H. 2010, A&A, 520, A43Ormel, C. W., Shi, J.-M., & Kuiper, R. 2015, MNRAS, 447, 3512 Owen, J. E. & Kollmeier, J. A. 2017, MNRAS, 467, 3379Paardekooper, S.-J., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401,1950Paardekooper, S.-J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293Paardekooper, S.-J. & Mellema, G. 2006, A&A, 459, L17Paardekooper, S.-J. & Mellema, G. 2008, A&A, 478, 245Paardekooper, S.-J. & Papaloizou, J. C. B. 2009, MNRAS, 394, 2283Petersen, M. R., Julien, K., & Stewart, G. R. 2007a, ApJ, 658, 1236Petersen, M. R., Stewart, G. R., & Julien, K. 2007b, ApJ, 658, 1252Pfeil, T. & Klahr, H. 2019, ApJ, 871, 150Piso, A.-M. A., Youdin, A. N., & Murray-Clay, R. A. 2015, ApJ, 800, 82Popovas, A., Nordlund, Å., & Ramsey, J. P. 2019, MNRAS, 482, L107Popovas, A., Nordlund, Å., Ramsey, J. P., & Ormel, C. W. 2018, MNRAS, 479,5136Raettig, N., Lyra, W., & Klahr, H. 2013, ApJ, 765, 115Rafikov, R. R. 2002, ApJ, 572, 566Rein, H. & Liu, S.-F. 2012, A&A, 537, A128Rein, H. & Spiegel, D. S. 2015, MNRAS, 446, 1424Ruden, S. P. & Pollack, J. B. 1991, ApJ, 375, 740Rüdiger, G., Arlt, R., & Shalybkov, D. 2002, A&A, 391, 781Semenov, D., Henning, T., Helling, C., Ilgner, M., & Sedlmayr, E. 2003, A&A,410, 611Stone, J. M. & Balbus, S. A. 1996, ApJ, 464, 364Stone, J. M. & Norman, M. L. 1992, ApJS, 80, 753Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257Tanaka, H. & Ward, W. R. 2004, ApJ, 602, 388Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47Ward, W. R. 1986, Icarus, 67, 164Ward, W. R. 1991, in Lunar and Planetary Inst. Technical Report, Vol. 22, Lunarand Planetary Science ConferenceWard, W. R. 1997, ApJ, 482, L211Yang, L. T. & Brent, R. P. 2002, in Proceedings of the Fifth Conference onAlgorithm and Architectures for Parallel Processing., ed. W. Zhou, X. Chi,A. Goscinski, & G. Li (IEEE), 324–328Yu, C., Li, H., Li, S., Lubow, S. H., & Lin, D. N. C. 2010, ApJ, 712, 198Zhu, Z., Hartmann, L., Nelson, R. P., & Gammie, C. F. 2012, ApJ, 746, 110
Article number, page 17 of 19 & A proofs: manuscript no. paper -5e-05-4e-05-3e-05-2e-05-1e-05 0 1e-05 2e-05 3e-05 4e-05 15 20 25 30 35 40 45 50 55 60 t o t a l t o r que Γ [ a J Ω J ] no. of orbits [ P orb ]Hot protoplanetCold protoplanet κ BL , instantaneous L κ BL , increasing L Fig. B.1.
Torque evolution in the κ BL -disk with instantaneously in-creased luminosity of the protoplanet L (solid black curve; same as inFig. 2) compared to the κ BL -disk with smoothly increasing L (dashedblue curve). Even before L reaches its maximum value in the latter case,the curves start to overlap. After t = P orb , the agreement is almost ex-act. Appendix A: Modifications of the numerical scheme
We implemented Eqs. (3) and (4) into fargo d using their dis-crete 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 sub-step. Using the same notation as in Bitsch et al. (2013), the re-lation between the temperature and radiative energy at t + ∆ t is T n + = η + η E n + , (A.1)and we redefine η = T n + ∆ t κ P c V σ ( T n ) + ∆ t Q (cid:15) − indep ρ c V + ∆ t κ P c V σ ( T n ) + ∆ t ( γ − ∇ · v , (A.2) η = ∆ t κ P c V c + ∆ t κ P c V σ ( T n ) + ∆ t ( γ − ∇ · v . (A.3)There are two changes with respect to Bitsch et al. (2013).First, the compressional heating is included via the last term inthe denominator of Eqs. (A.2) and (A.3). Second, the Q (cid:15) − indep term is a sum of all heat sources that do not depend on (cid:15) (or T );in our case Q (cid:15) − indep = Q visc + Q art + Q acc . We note that whennecessary, the stellar irradiation term can also be easily includedin Q (cid:15) − indep but it is neglected in this work. Appendix B: Supporting simulations
In this appendix, we summarize several additional simulationsdesigned to confirm the robustness of our conclusions.
Appendix B.1: Simulation with a smoothly increasingluminosity of the protoplanet
In our main simulations, we usually start the phase with accre-tion heating abruptly, by instantaneously increasing the luminos-ity of the protoplanet from L = τ =
100 kyr. Such a sudden appearanceof a strong heat source might produce an unexpected behaviourand instabilities by itself. In order to exclude any undesirable be-haviour, we repeat the κ BL -disk simulation with accretion heating -6e-05-5e-05-4e-05-3e-05-2e-05-1e-05 0 1e-05 2e-05 3e-05 4e-05 6 8 10 12 14 t o t a l t o r que Γ [ a J Ω J ] no. of orbits [ P orb ]Hot protoplanetCold protoplanet κ BL , higher resol. Fig. B.2.
Torque evolution obtained in the κ BL -disk simulation withan increased azimuthal resolution of 2764 cells. The time span of theindividual phases is shortened to save computing time. The instabilityin the presence of the accretion heating is, however, recovered again. of the protoplanet, but now we linearly increase L from zero at t = P orb to its maximal value over the time interval of 10 P orb .Fig. B.1 shows the measured torque evolution. Clearly, thegradual increase of L has no impact on the final character ofthe torque, the oscillations related to the flow reconfigurationsinevitably appear. Appendix B.2: Simulation with an increased azimuthalresolution
The resolution in our simulations is motivated by works of Legaet al. (2014) and Eklund & Masset (2017). Although we do notperform extended convergence tests of our own, the resolutionshould be su ffi cient to recover a realistic torque value and alsoa realistic advection-di ff usion redistribution of the hot gas nearthe protoplanet.However, once the circumplanetary flow becomes unstable,it is no longer clear if the chosen resolution is su ffi cient. Forexample, one might argue that the coverage of the Hill sphereby the grid cells in the azimuthal direction is too poor. Here wepresent an experiment in which we double the number of thegrid cells in the azimuthal direction in order to have the samecoverage of the Hill sphere in all directions. We perform the κ BL -disk simulation again, however, we shorten the phase withoutaccretion heating to 10 P orb and the phase with accretion heatingto 5 P orb .The result is shown in Fig. B.2 and demonstrates that thetoque oscillations are recovered even when the increased reso-lution is used. However, the amplitude of the torque oscillationsslightly changes, implying that the resolution dependence shouldbe explored more carefully in the future. Appendix B.3: Simulation with the unmodified Bell & Linopacity table
The simulations of the κ BL -disk presented in this paper are per-formed with a simplified opacity law of Bell & Lin (1994) (ex-plained in Sect. 2.2). Here we test how the results change if theunmodified opacity law κ fullBL is used and the dependence on thelocal values of T and ρ is retained.Fig. B.3 compares the torque evolution in our standard κ BL -disk with that in a κ fullBL -disk. Clearly, the instability occurs inboth disks, regardless of whether or not the input values for theopacity function are azimuthally averaged. The only di ff erenceis in the torque amplitude which is larger in the κ fullBL -disk. Article number, page 18 of 19. Chrenko and M. Lambrechts: Oscillatory migration of accreting protoplanets driven by a 3D distortion of the gas flow -8e-05-6e-05-4e-05-2e-05 0 2e-05 4e-05 6e-05 15 20 25 30 35 40 45 50 55 60 t o t a l t o r que Γ [ a J Ω J ] no. of orbits [ P orb ]Hot protoplanetCold protoplanet κ BLfull κ BL Fig. B.3.
Comparison of the torque evolutions obtained in the κ BL -disk(solid black curve; same as in Fig. 2) and the κ fullBL -disk (dashed bluecurve). The first model inputs azimuthally-averaged values of ρ and T to the opacity function of Bell & Lin (1994), while the latter uses localvalues of ρ and T . Both cases lead to the instability of the circumplane-tary flow, but they di ff er in the amplitude of the torque oscillations. -8e-05-6e-05-4e-05-2e-05 0 2e-05 4e-05 6e-05 20 25 30 35 40 t o t a l t o r que Γ [ a J Ω J ] no. of orbits [ P orb ]Hot protoplanetCold protoplanetFargoca, κ BLfull
Fargo3D, κ BLfull
Fig. B.4.
Comparison of the torque evolutions obtained with two in-dependent codes fargo d (solid black curve) and fargoca (dashed bluecurve). The unmodified κ fullBL opacity table of Bell & Lin (1994) was usedin this case. The increased amplitude occurs because if T locally rises,so does the material opacity ( κ fullBL ∝ T in the given disk re-gion). Subsequently, the radiative cooling of the hot gas becomesless e ffi cient and T rises even more. As a result, the underdenseperturbations related to any temperature excess are more pro-nounced, leading to the larger amplitude of the torque oscilla-tions. Appendix B.4: Code comparison
With the aim to confirm that our implementation of the energyequations in fargo d is correct and also that the instability ofthe circumplanetary flow does not arise due to numerical arte-facts, we present a comparison simulation obtained with an in-dependent and well-tested radiation hydrodynamic code fargoca (Lega et al. 2014). The simulation is performed with the unmod-ified κ fullBL opacity law of Bell & Lin (1994).Fig. B.4 compares the torque evolution found using our codewith the one obtained with fargoca . We can see that the con-verged torque for the cold protoplanet is in a satisfactory agree-ment.When the accretion heating is initiated, the same oscillatorytrend is observed with both codes. The curves overlap at first;later they start to depart in terms of the oscillation phase. How-ever, the amplitude remains the same.Since the converged torques are in agreement and the in-stability is recovered, we conclude that the di ff erences that we identify in Fig. B.4 arise only because our radiation module in fargo d relies on a slightly di ff erent numerical scheme (see Ap-pendix A) compared to fargoca . Appendix C: Vorticity equation in the corotatingframe
For reader’s convenience, we provide a step-by-step derivationof the vorticity equation. Starting with Eq. (2), we apply thecurl on both sides. In the following, we neglect the viscous term(large Reynolds number limit).The following identities of the vector calculus will beutilised:( a · ∇ ) a = ∇ ( a · a ) + ( ∇ × a ) × a , (C.1) ∇ × ( a × b ) = a ( ∇ · b ) − b ( ∇ · a ) + ( b · ∇ ) a − ( a · ∇ ) b , (C.2) ∇ · ( ∇ × a ) = , (C.3) ∇ × ( ∇ f ) = , (C.4)where the last identity holds for scalar functions that are at leasttwice continuously di ff erentiable.The curl of the advection term yields ∇ × [( v · ∇ ) v ] = ∇ × (cid:34) ∇ ( v · v ) + ( ∇ × v ) × v (cid:35) = ∇ × [ ω r × v ] , (C.5)where we used Eq. (C.1) in writing the first equality, Eq. (C.4)to remove the ∼ v · v term and we defined the relative vorticity inthe corotating frame ω r = ∇ × v . Using Eqs. (C.2) and (C.3), wefurther obtain ∇ × [( v · ∇ ) v ] = ω r ( ∇ · v ) + ( v · ∇ ) ω r − ( ω r · ∇ ) v , (C.6)The curl of the pressure term leads to ∇ × (cid:32) ∇ P ρ (cid:33) = − ρ ∇ ρ × ∇ P , (C.7)because ∇ × ( ∇ P ) = (Eq. C.4).When dealing with the gravitational term, it is useful torealise that the centrifugal acceleration can be expressed as Ω × ( Ω × r ) = ∇ Φ c , with the centrifugal potential Φ c = − r ⊥ Ω .The curl of a combined force term, ∇ × [ ∇ ( Φ + Φ c )], is zero ow-ing to Eq. (C.4).Finally, we take the curl of the Coriolis acceleration: ∇ × (2 Ω × v ) = Ω ( ∇ · v ) − ( Ω · ∇ ) v ] , (C.8)where we removed terms ∼∇ · Ω , ∼∇ Ω because Ω in our simu-lations is constant.Recollecting all the terms, we can write the relative vorticityequation ∂ ω r ∂ t + ( v · ∇ ) ω r = D ω r D t = ( ω a · ∇ ) v − ω a ( ∇ · v ) + ∇ ρ × ∇ P ρ , (C.9)where we defined the absolute vorticity in the inertial frame ω a = ω r + Ω ..