# Planetary Torque in 3D Isentropic Disks

Jeffrey Fung, Frédéric Masset, Elena Lega, David Velasco

► To cite this version:

Jeffrey Fung, Frédéric Masset, Elena Lega, David Velasco. Planetary Torque in 3D Isentropic Disks. The Astronomical Journal, 2017, 153 (3), pp.124. 10.3847/1538-3881/153/3/124. hal-03252749

HAL Id: hal-03252749

https://hal.science/hal-03252749v1

Submitted on 25 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Planetary Torque in 3D Isentropic Disks

Jeffrey Fung<fnref n="1" />,<fnref n="4" />, Frédéric Masset<fnref n="2" />, Elena Lega<fnref n="3" />, and David Velasco<fnref n="2" />

Received 2016 October 28; revised 2016 December 19; accepted 2017 January 10; published 2017 February 21

## Abstract

Planetary migration is inherently a three-dimensional (3D) problem, because Earth-size planetary cores are deeply embedded in protoplanetary disks. Simulations of these 3D disks remain challenging due to the steep resolution requirements. Using two different hydrodynamics codes, FARGO3D and PEnGUIn, we simulate disk–planet interaction for a one to five Earth-mass planet embedded in an isentropic disk. We measure the torque on the planet and ensure that the measurements are converged both in resolution and between the two codes. We find that the torque is independent of the smoothing length of the planet’s potential ($r_s$), and that it has a weak dependence on the adiabatic index of the gaseous disk ($\gamma$). The torque values correspond to an inward migration rate qualitatively similar to previous linear calculations. We perform additional simulations with explicit radiative transfer using FARGOCA, and again find agreement between 3D simulations and existing torque formulae. We also present the flow pattern around the planets that show active flow is present within the planet’s Hill sphere, and meridional vortices are shed downstream. The vertical flow speed near the planet is faster for a smaller $r_s$ or $\gamma$, up to supersonic speeds for the smallest $r_s$ and $\gamma$ in our study.

**Key words:** accretion, accretion disks – circumstellar matter – methods: numerical – planet–disk interactions – planets and satellites: formation – protoplanetary disks

## 1. Introduction

Newly born planets interact gravitationally with their natal circumstellar disks. As the planet’s tidal force exerts a torque on the disk, the back-reaction from the disk also torques the planet, causing it to migrate. The study of planet migration is one of the links that connects the initial formation of planets to their final positions in their planetary systems, and is therefore essential to explaining the statistical distribution of planets (e.g., Petigura et al. 2013; Dressing & Charbonneau 2015; Silburt et al. 2015). Calculations of planet migration rates have been done extensively for planets interacting with razor-thin, two-dimensional (2D) disks (for a review, see Baruteau & Masset 2013 and references therein), and, to a lesser extent, with more realistic, three-dimensional disks (3D; e.g., Tanaka et al. 2002; Bate et al. 2003; D’Angelo et al. 2003; Masset et al. 2006; D’Angelo & Lubow 2010; Uribe et al. 2011; Lega et al. 2014; Benítez-Llambay et al. 2015).

Recently, Fung et al. (2015, hereafter FAW15) reported a torque on an Earth-size planet embedded in an isothermal disk that significantly deviates from linear estimates, such as those by Tanaka et al. (2002). This opens the possibility that the co-orbital region may behave in a way previously unanticipated, thus generating a new component in the torque. Their results, however, also show insufficient resolution near the planet, leading to the need to exclude the planet’s vicinity for their torque calculation. The existence of this new torque component therefore demands verification.

Also demonstrated by recent results in the literature, including those of Benítez-Llambay et al. (2015) and FAW15, is that in 3D disks, the density structure close to the planet, within the scale of the planet’s Hill radius, can potentially have an overwhelming influence on the planetary torque. It is therefore essential that the thermal properties of the

gas in this area be properly modeled. Protoplanetary disks models commonly assume a locally isothermal structure. Although this may be appropriate for the global disk, it is unsuitable for the gas surrounding the planet. In reality, the planet’s atmosphere is opaque to its own thermal emission and cools less efficiently than the disk, and so it should more resemble a heated envelope within which convection is expected to operate, generating a constant entropy profile. This implies that an isentropic equation of state may be more applicable to the density structure close to the planet. Additionally, the hydrostatic profile of an isothermal gas inside the planet’s gravitational potential is exponential, which makes achieving numerical convergence challenging. As we will show later in this paper, convergence is more readily achieved when the adiabatic index, $\gamma$, is larger.

This paper aims to obtain a converged measurement of the torque on an Earth-size planet embedded in a 3D isentropic disk. We will investigate possible dependencies on both the adiabatic index, $\gamma$, and the smoothing length of the planet’s gravitational potential, $r_s$. To identify potential code bias, we simulate identical models with two hydrodynamics codes: FARGO3D and PEnGUIn. To further identify possible discrepancies between 2D and 3D disks, we also present additional simulations with explicit radiative transfer using FARGOCA.

Section 2 contains code descriptions and simulation parameters. Section 3 presents our results for both the torque and the flow structure around the planet. Section 4 concludes and discusses implications of our results.

## 2. Numerical Method

The simulations are performed in spherical coordinates, where $r$, $\phi$, and $\theta$ denote the usual radial, azimuthal, and polar coordinates. For convenience, we also denote $R \equiv r \sin\theta$ and $z \equiv r \cos\theta$ as the cylindrical radial and vertical coordinates.

<footnote n="1">¹Department of Astronomy, University of California at Berkeley, Campbell Hall, Berkeley, CA 94720-3411, USA; jeffrey.fung@berkeley.edu</footnote>

<footnote n="2">²Instituto de Ciencias Físicas, Universidad Nacional Autónoma de México, Av. Universidad s/n, 62210 Cuernavaca, Mor., Mexico</footnote>

<footnote n="3">³Université de la Côte d'Azur, Observatoire de la Côte d'Azur, CNRS, Laboratoire Lagrange UMR 7293, Nice, France</footnote>

<footnote n="4">⁴ NASA Sagan Fellow.</footnote>
FARGO3D and PEnGUIn solve the compressible Euler equations:

$$ \frac{D\rho}{Dt} = -\rho(\nabla \cdot \mathbf{v}), \qquad (1) $$

$$ \frac{D\mathbf{v}}{Dt} = -\frac{1}{\rho}\nabla p - \nabla\Phi, \qquad (2) $$

where $D/Dt$ is the Lagrangian derivative, $\rho$ is the density, $p$ the gas pressure, $\mathbf{v}$ the velocity, and $\Phi$ the gravitational potential of the central star and the planet, representing the only external force we consider in this study.

For our model with explicit radiative transfer, we include the following energy equation:

$$ \frac{De}{Dt} = -e(\nabla \cdot \mathbf{v}) - p(\nabla \cdot \mathbf{v}) - \nabla \cdot \mathbf{F} + Q^{+}, \qquad (3) $$

where $e$ is the thermal energy density: $e = \rho c_v T$, $c_v$ the specific heat at constant volume and $T$ the gas temperature. The term $p(\nabla \cdot \mathbf{v})$ represents compressional heating, $Q^+ = (\nabla \nabla) \cdot \mathbf{v}$ viscous heating (with $\nabla$ the viscous stress tensor, as in Mihalas & Weibel Mihalas 1984) and $\mathbf{F}$ is the radiative flux. We have chosen a disk with a constant viscosity<fnref n="5" /> $\nu = 10^{-5} R_p^2 \Omega_p$, where $R_p$ is the orbital radius of the planet's orbit, and $\Omega_p$ the planet's orbital frequency. The radiative flux is calculated in the framework of the flux-limited diffusion approximation (Levermore & Pomraning 1981),

$$ \mathbf{F} = -\frac{c_f 4 a_R T^3}{\rho \kappa} \nabla T \qquad (4) $$

where $f_\lambda$ is a flux-limiter (Kley 1989), $c$ is the speed of light, $a_R$ the radiation constant, and $\kappa$ the Rosseland mean opacity, set to be $0.01 \text{ cm}^2 \text{ g}^{-1}$ (or $3.28 \times 10^3$ in code units).

We note that, for our isentropic models where viscous heating is not a concern (unlike the radiative case), we have deliberately chosen not to include viscosity because we consider it a more physical choice. Rather than being viscous, real protoplanetary disks are turbulent, and it is unclear whether this turbulence translates to viscosity when considering detail flow structures on the scale of the planet's Hill radius. In fact, viscosity could easily generate unrealistic behavior, such as preventing circumplanetary flow by removing angular momentum too rapidly.

The simulations are performed in a frame centered on the star, and we fix the planet on a circular orbit in the disk midplane. Therefore, $\Phi$ is:

$$ \Phi = -\frac{GM}{1+q} \left[ \frac{1}{r} + \frac{q}{\sqrt{r^2 + R_p^2 - 2RR_p \cos \phi' + r_s^2}} - \frac{qR \cos \phi'}{R_p^2} \right], \qquad (5) $$

where $G$ is the gravitational constant, $M = M_\ast + M_p$ the total mass of the star and the planet, $q = M_p/M_\ast$ the planet-to-star mass ratio, $r_s$ the smoothing length of the planet's potential, and $\phi' = \phi - \phi_p$ denotes the azimuthal separation from the

planet. We set $GM = 1$ and $R_p = 1$, so that the Keplerian velocity and frequency $v_k = \sqrt{GM/r}$ and $\Omega_k = \sqrt{GM/r^3}$ both equal 1 at the planet's orbit. We also label the planet's orbital speed as $v_p$ for convenience. The third term in the bracket is the indirect potential due to the non-inertial frame.

For comparison with FAW15, we choose $q = 1.5 \times 10^{-5}$, or about five Earth masses for a Solar-mass star, for most of our simulations, as well as an additional few with $q = 3 \times 10^{-6}$, or about one Earth mass. The Hill radius of a planet is:

$$ r_H = R_p \left(\frac{q}{3}\right)^{\frac{1}{3}}, \qquad (6) $$

which comes to $r_H \approx 0.017 R_p$ for a five Earth-mass ($M_\oplus$) planet. We set $r_s$ to be a small fraction of $r_H$, and look into torque's dependence on it. Our choices for $r_s$ are between 3% and 10% of $r_H$. Realistically, $r_s$ should correspond to the physical size of the planet, which, for this example, is about 0.4% of $r_H$ for the Earth. However, it is not numerically feasible to use such a small $r_s$, as the resolution required would be prohibitively high.

Completing our set of equations is the isentropic equation of state:

$$ p = \frac{c_0^2 \rho_0}{\gamma} \left(\frac{\rho}{\rho_0}\right)^{\gamma}, \qquad (7) $$

where $c_0$ is the adiabatic sound speed when $\rho = \rho_0$, and the normalization $\rho_0$ is set to 1. The sound speed is calculated from $c_s = \sqrt{\gamma p/\rho}$. In all of our simulations, we fix $c_0 = 0.03$ and scale our density profile such that $\rho = \rho_0$ at the planet's location. As a result, the sound speed near the planet is always 3% of the local Keplerian speed, regardless of $\gamma$. This is the same sound speed as the one used by FAW15. We simulate disks with three different values of $\gamma$: 1.2, 1.4, and 1.67.

For the model with radiative transfer, we instead use the ideal gas equation of state:

$$ p = \frac{R_{\text{gas}} \rho T}{\mu} \qquad (8) $$

with mean molecular weight $\mu$, $\mu = 2.3$, and gas constant $R_{\text{gas}}$. The sound speed is calculated using $\gamma = 1.4$. The dependence of the torque on $\gamma$ for radiative models was previously investigated by Bitsch et al. (2013).

## 2.1. Initial and Boundary Conditions

The disks are initialized assuming hydrostatic equilibrium. The density profile is:

$$ \rho = \rho_0 \left[ \left( \frac{R}{R_p} \right)^{-(\beta+\frac{3}{2})\frac{2(\gamma-1)}{\gamma+1}} - \frac{GM(\gamma-1)}{c_0^2} \left( \frac{1}{R} - \frac{1}{r} \right) \right]^{\frac{1}{\gamma-1}}, \qquad (9) $$

where $\beta$ defines the surface density profile $\Sigma \propto R^{-\beta}$. In this work, we choose $\beta = 3/2$. This produces a constant initial vortensity profile, $(\nabla \times \mathbf{v}_k)/\Sigma$, which, in 2D disks, would imply a net zero horseshoe drag (Ward 1991). Masset & Benítez-Llambay (2016) have shown that the same applies in 3D, as long as vorticity is kept in the vertical direction (i.e., no vortex tilting). With both a constant entropy and vortensity,

<footnote n="5">⁵ For the radiative disk, the rhs of Equation (2) has an additional term corresponding to the divergence of the viscous stress tensor $(-\frac{1}{\rho}\nabla \cdot \mathbf{T})$.</footnote>
none of the known sources of corotation torque are present in our model.

The orbital frequency of the disk in hydrostatic equilibrium is modified by the radial pressure gradient:

$$ \Omega = \sqrt{\Omega_k^2 + \frac{1}{r\rho} \frac{\partial p}{\partial r}}, \quad (10) $$

whereas the radial and polar velocities are zero. The planet is introduced to the disk gradually, where its mass increases to the desired value over the first orbit.

Our simulation domain spans $0.7R_p$ to $1.3R_p$ in the radial, and the full $2\pi$ in azimuth. Because our density profile is polytropic, it falls to negative values when

$$ z > z_{\text{lim}} \equiv \sqrt{\frac{2}{\gamma - 1} \frac{c_0}{\Omega_k} R^{-\left(\beta + \frac{3}{2}\right) \frac{\gamma - 1}{\gamma + 1}}.} \quad (11) $$

Therefore, we determine the top polar boundary by setting it below $arctan(z_{\text{lim}}/R)$ for all $R$ within our simulation domain, which is 0.087 radian above the midplane when $\gamma = 1.2$, 0.066 when $\gamma = 1.4$, and 0.048 when $\gamma = 1.67$. For the radiative simulations, the top polar boundary is placed at 0.1 radian above the midplane. In all models, we set our bottom boundary at the disk midplane.

For our radial boundaries, we apply a fixed boundary condition where the disk variables remain at their initial values. An additional wave-killing zone (de Val-Borro et al. 2006) is implemented in the radiative simulations to prevent the reflection of density waves. The polar boundaries are reflective, both at the top, to prevent mass from entering or leaving the domain, and at the bottom, to ensure the disk is symmetric across the midplane.

## 2.2. Code Descriptions and Simulation Resolution

We use three different resolutions to test the numerical convergence of our results; a “low” resolution where the planet’s Hill radius $r_H$ is resolved by ~30 cells, a “medium” resolution where it is doubled to ~60 cells, and a “high” resolution with resolution redoubled to ~120 cells. Our codes use different methods to achieve these resolutions in the vicinity of the planet: nested grid in FARGO3D, and non-uniform grid geometry in PEnGUIn and FARGOCA.

In our three codes, the reference frame co-rotates with the planet, but the Coriolis force is not computed as an explicit source term, as shown by Equation (2). Rather, it is absorbed into the conservative form of the angular momentum equation (Kley 1998).

In the following, we give a brief description of these codes.

### 2.2.1. FARGO3D

The code FARGO3D (Benítez-Llambay & Masset 2016) is here used with a newly implemented nested grid capability. On top of a base grid that has the boundaries specified in Section 2.1, we use a hierarchy of nested grids with a doubling-up of the resolution between successive grid levels. The limits of our nested grids are given in Table 1. Within each grid, the hydrodynamical solver described by Benítez-Llambay & Masset (2016) is used, with orbital advection deactivated. Boundary conditions of the nested meshes are imposed by performing a trilinear interpolation of the underlying coarser mesh in a three-cell wide layer of ghost zones surrounding the

<table><caption>**Table 1**
Position of the Nested Grids</caption><thead><tr><th>φ<sub>min</sub></th><th>r<sub>min</sub>/R<sub>p</sub></th><th>θ<sub>min</sub></th><th>φ<sub>max</sub></th><th>r<sub>max</sub>/R<sub>p</sub></th><th>θ<sub>max</sub></th><th>level</th></tr></thead><tbody><tr><td>-1.0</td><td>0.89</td><td>1.54</td><td>1.0</td><td>1.11</td><td>π/2</td><td>1</td></tr><tr><td>-0.2</td><td>0.91</td><td>1.545</td><td>0.2</td><td>1.09</td><td>π/2</td><td>2</td></tr><tr><td>-0.1</td><td>0.95</td><td>1.553</td><td>0.1</td><td>1.05</td><td>π/2</td><td>3</td></tr><tr><td>-0.04</td><td>0.97</td><td>1.556</td><td>0.04</td><td>1.03</td><td>π/2</td><td>4</td></tr><tr><td>-0.025</td><td>0.98</td><td>1.559</td><td>0.025</td><td>1.02</td><td>π/2</td><td>5</td></tr><tr><td>-0.018</td><td>0.985</td><td>1.562</td><td>0.018</td><td>1.015</td><td>π/2</td><td>6</td></tr></tbody></table>

**Note.** Levels 1 to 4 are used in the “low” resolution runs, 1 to 5 in the “medium” resolution runs, and 1 to 6 in the “high” resolution runs. This table lists the intended limits, set by the user. The actual limits differ slightly, as they have to fall on the interfaces between cells of the coarser level.

mesh, except for the boundaries at the midplane ($\theta = \pi/2$), which are reflective. Upon integration on a given level (hereafter, fine level), the information is communicated to the underlying coarser level (hereafter, coarse level) in two ways: (i) in the outermost two-cell wide contour of the coarse level covered by the fine one, the different hydrodynamics quantities are replaced by their averaged values from the fine level, and (ii) the fluxes of mass and momentum on the contour of the fine level are used to update the quantities on the coarse level, which ensures that the code conserves mass and angular momentum to machine accuracy. The integration is done recursively across the whole hierarchy of nested grids. If a given level is advanced in time with a time step $\Delta t$, the next finer level is either also advanced over a time step $\Delta t$, or it is advanced twice with a time step $\Delta t/2$. Which of these two possibilities is chosen depends on which one yields the largest advance in time for a given computational cost. Our implementation runs on graphical processing units (GPUs) and is parallelized using a message passing interface (MPI). The details of our implementations will be presented elsewhere.

### 2.2.2. PEnGUIn

PEnGUIn (Piecewise Parabolic Hydro-code Enhanced with GPU Implementation) is a Lagrangian, dimensionally split, shock-capturing hydrodynamics code that runs on GPUs (Fung 2015).

Different from FARGO3D, PEnGUIn achieves high resolution using a non-uniform grid geometry. We assign small cell sizes near the planet, and continuously increase the cell size as the distance from the planet grows. The prescription is as follows:

$$ \Delta = \frac{x_{\text{max}}^a}{Na} x^{1-a} + \Delta_{\text{min}}, \quad (12) $$

$$ a = \frac{x_{\text{max}}}{N \Delta_{\text{max}}}, \quad (13) $$

where $\Delta$ is the cell size, $x$ is the distance away from the planet, $x_{\text{max}}$ is the maximum value of $x$ (i.e., the domain boundary), $N$ is the number of cells between 0 and $x_{\text{max}}$, and $\Delta_{\text{min}}$ and $\Delta_{\text{max}}$ are the assigned cell sizes at $x = 0$ and $x = x_{\text{max}}$ respectively. For our low-, medium-, and high-resolution runs, we set $\Delta_{\text{min}} = 5.0 \times 10^{-4}$, $2.5 \times 10^{-4}$, and $1.25 \times 10^{-4}$ respectively. We always fix $\Delta_{\text{max}}$ at 0.002 in the radial and polar directions, and 0.015 in the azimuthal direction.
### 2.2.3. FARGOCA

The code FARGO with Colatitude Added (FARGOCA) (Lega et al. 2014) is based on the FARGO code (Masset 2000) extended to three-dimensions with the additional introduction of an energy equation to provide a realistic modeling of radiative effects (Kley et al. 2009). The fluid equations are solved using finite differences with a time-explicit-implicit multistep procedure. Precisely, concerning the energy equation, we first update the energy by explicit integration of the compressional heating term, and in a separate step we integrate the viscous heating and the radiative diffusion terms. In this second step, we follow the backward Euler method, which is an unconditionally stable implicit method solved with a standard SOR (Successive Overrelaxation Reduction) solver. The code is parallelized using a hybridization of MPI between the nodes and OpenMP on shared memory multi-core processors. High resolution is achieved, as in PEnGUIn, using the nonuniform grid geometry with the prescription detailed above. Because hydrodynamical 3D calculations of radiative disks are very expensive in computational time, we have moderately lower resolution with respect to the values used for isentropic disks with both FARGO3D and PEnGUIn. We set in the following for our low-, medium-, and high-resolution radiative runs: $\Delta_{min} = 1.2 \times 10^{-3}$, $6 \times 10^{-4}$, and $3 \times 10^{-4}$, respectively.

## 3. Results

### 3.1. Planetary Torque

We run each simulation to five orbits, at which point a steady torque on the planet has been established, but the planet has yet to create nonlinear modifications to the disk structure, such as planetary gaps. We defer a discussion on the long-term effects to Section 4. We measure the net torque on the planet without excluding any region, and label the torque $\Gamma$, normalizing it by $\Gamma_0 \equiv \Sigma_0 R_p^4 \Omega_p^2 q^2 (c_0/v_p)^{-2}$ (this also applies to the radiative model, where the equilibrium temperature at the planet's position has the same sound speed as $c_0$).

Convergence requires higher resolution for smaller $r_s$ or $\gamma$ values. For instance, we find that our low-resolution setup is typically sufficient when $r_s = 0.1 r_H$, but high resolution is needed when $r_s < 0.05 r_H$. Typically, $r_s$ needs to be resolved by at least three cells when $\gamma \gtrsim 1.4$, and about four to five cells when $\gamma = 1.2$. We achieved convergence to within a few percent for nearly all of our models, both in resolution for FARGO3D and PEnGUIn individually, and in code comparison between them. Figure 1 plots the torque of our $\{r_s/r_H, \gamma\} = \{0.03, 1.2\}$ model, as a function of time. This model has the lowest $r_s$ and $\gamma$ values in our parameter space, and shows the highest level of fluctuation; nonetheless, the agreement between PEnGUIn and FARGO3D is within 5%. In the rest of this paper, we will only present the converged torque measurements, which are either from medium- or high-resolution simulations, with the medium-resolution results verified by shorter high-resolution runs.

In Figure 2, we find that the torque has no dependence on $r_s$ for the range of values we considered ($r_s = 0.03-0.1 r_H$). The models shown are of $\gamma = 1.4$, but the same is found for other $\gamma$ values as well. In Figure 3, we plot models with different $\gamma$s, and observe a weak trend where the torque is more negative for a smaller $\gamma$. In PEnGUIn simulations, we measure $-2.5, -2.3$, and $-2.2 \Gamma_0$ for $\gamma = 1.2, 1.4$, and $1.67$, respectively; FARGO3D simulations show the same trend with torques of

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Torque on a $5 M_\oplus$ planet as a function of time for different resolutions, with both FARGO3D and PEnGUIn. All models shown have $\{r_s/r_H, \gamma\} = \{0.03, 1.2\}$. Solid, dashed, and dashed-dotted-dotted lines correspond to “high,” “medium,” and “low” resolutions, as defined in Section 2.2. FARGO3D simulations are in blue, cyan, and green, whereas PEnGUIn ones are in red, magenta, and orange. Data points are time-averaged over 0.1 orbit. At high resolution, both codes converge to the same torque value after five orbits, $\sim -2.5 \Gamma_0$, and the agreement between them is within 5%.</figcaption></figure>

$-2.6, -2.4$, and $-2.2$, respectively. Figure 6 plots the torque density $d\Gamma/dM$ (torque per mass at each annulus) of these models, with the 3D isothermal model of D’Angelo & Lubow (2010) overlayed (interpolated between their $\Sigma \propto r^{-1}$ and $r^{-2}$ models). We find a torque arising in the planet’s co-orbital region as $\gamma$ decreases, which largely accounts for the observed trend in $\gamma$.

This corotation torque is not predicted by the linear calculations of Tanaka et al. (2002). This is likely because the $5 M_\oplus$ planet we choose is about half the disk’s thermal mass:

$$M_{\text{thermal}} = \left(\frac{c_0}{v_p}\right)^3 M_*, \quad (14)$$

and $M_p \sim 0.56 M_{\text{thermal}}$. It is known that, in this regime, the flow pattern around the planet begins to deviate from linear calculations (e.g., Korycansky & Papaloizou 1996), which may well lead to modifications to the linear torque. To confirm this, we perform the same simulations, but for a $1 M_\oplus$ planet, which is highly sub-thermal ($M_p \ll M_{\text{thermal}}$). The results are shown in Figure 4. For this smaller planet, the trend in $\gamma$ is reversed, and becomes compatible with the isothermal linear estimate. Figure 7 plots the torque density for some of these models, clearly showing the absence of any co-orbital component. The model by D’Angelo & Lubow (2010) is also in the highly sub-thermal regime, and we clearly see that their torque density follows the trend in $\gamma$ much better in Figure 7 than in Figure 6. Therefore, we conclude that the corotation torque in the $5 M_\oplus$ models is nonlinear in nature.

Corotation torques in inviscid disks are expected to saturate over a few libration times. This cannot be seen with our five-orbit run time, and much longer simulations are beyond the capacity of our computational resources. We do, however, note that torque saturation is a prediction for 2D linear flow, and it is unclear whether the same applies for 3D nonlinear flow.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Torque on a 5 $M_\odot$ planet as a function of time for different values of the smoothing length $r_s$. The left panel plots simulation results from FARGO3D, and right panel from PEnGUIn. All models shown have $\gamma = 1.4$. Overall, we find the torque to be independent of $r_s$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Torque on a 5 $M_\odot$ planet as a function of time for different values of the adiabatic index $\gamma$. All models shown have $r_s = 0.07 r_\text{H}$. Both codes show that the torque has a weak dependence on $\gamma$.</figcaption></figure>

Although our two codes are generally in excellent agreement, we do find torque measurements with FARGO3D to have stronger fluctuations than PEnGUIn. This is in line with the more turbulent flow pattern it finds, which we discuss in Section 3.2. We suspect that they are numerical artifacts generated from the interfaces between different mesh refinement levels. The level of fluctuation is sufficiently small that it does not affect our main conclusions.

Our torque measurements significantly deviate from the one reported by FAW15, which was -0.8 $\Gamma_0$. In terms of our parameters, their simulation uses $\{r_s/r_\text{H}, \gamma\} = \{0.1, 1\}$, and a resolution that translates to about three cells per $r_s$. Using our resolution study as a guideline, it is likely that their resolution was insufficient for a converged torque measurement. We expect, from extrapolating our results, that the converged torque measurement in their disk model should instead be less than -2.5 $\Gamma_0$. We have attempted to verify this prediction using our high-resolution grid with the same parameters $\{r_s/r_\text{H}, \gamma\} = \{0.1, 1\}$, but the resulting torque

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Similar to Figure 3, but for 1 $M_\odot$ planets. All models shown have $r_s = 0.07 r_\text{H}$. The dashed line indicates the linear torque computed by Tanaka et al. (2002) for an isothermal disk. Comparing to Figure 3, we find the trend in $\gamma$ is reversed, and is compatible with linear estimates in the isothermal limit.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Torque on a 5 $M_\odot$ planet as a function of time, for the radiative models from FARGOCA. Left panel: torque measured for different resolutions values. The “low,” “medium,” and “high” resolutions for the radiative case are defined in Section 2.2. The medium- and high-resolution runs are restarted from the low-resolution run after the torque has come to a steady state (about 20 orbits). We measure a final torque of ~0.55 $\Gamma_0$, in agreement with the expected analytic value. Right panel: torque measured for two different smoothing lengths $r_s$. Similar to Figure 2, we observe no dependence on $r_s$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Torque density from the same models as those in Figure 3, obtained at the end of our simulations. Also shown, as the dashed black curves, are the isothermal model of D’Angelo & Lubow (2010). The insets zoom in on the vicinity around the planet, revealing a small corotation torque that increases with decreasing $\gamma$. The two codes agree well, with the FARGO3D result showing more noise coming from the lower-resolution regions in the nested mesh.</figcaption></figure>

strongly fluctuates between -10 and $5 \Gamma_0$,<fnref n="6" /> which suggests that our resolution is still insufficient. Therefore, we caution the reader that 3D inviscid isothermal simulations may have a resolution requirement substantially more severe than commonly expected.

For the radiative simulations, more care has to be taken because the disk requires more time to adjust thermally. To obtain the final torque, we follow a three-step procedure: first, we perform a 2D ($r-\theta$), axisymmetric simulation of the disk (without a planet) to bring the disk to thermal equilibrium; second, expand the 2D disk along the azimuth into our 3D low resolution grid, introduce the planet, and continue the run for 50 orbits; finally, we interpolate the grid into medium or high resolution and continue for five additional orbits. Figure 5 plots

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Torque density from our 1 $M_\odot$, $r_s = 0.07 r_H$ PEnGUIn simulations. Similar to Figure 6, we also overlay the isothermal model of D’Angelo & Lubow (2010) here, as the dashed black curve. We find that the 1 $M_\odot$ models do not show any corotation torque, and the Lindblad torque is stronger for a larger $\gamma$.</figcaption></figure>

<footnote n="6">⁶ FAW15’s measurement and ours differ in two ways: FAW15 excluded a sphere of $0.5 r_H$ radius around the planet, whereas we made no excision; and FAW15’s simulation included a low level of viscosity, whereas ours is inviscid. Both of these differences stem from FAW15’s attempt to reduce numerical noise, which we avoid in order to obtain true numerically converged values.</footnote>
<figure><img src="image_10.png" /></figure>

**Figure 8.** Vertical velocity $v_z$ at an azimuthal slice across the planet's position. The planet is located at $\{R, z\} = \{R_p, 0\}$. Red (positive values) indicates velocities upward (away from the midplane), and blue (negative values) indicates downward (toward the midplane). $v_z$ is normalized to the initial disk sound speed at the planet's location $c_0$, and the color stretch scales with $\sqrt{v_z}$ to emphasize finer details. All panels show the $\{r_s/r_H, \gamma\} = \{0.03, 1.2\}$ model at high resolution for 5 $M_\oplus$ planets. The left panel is from FARGO3D at one orbit, the middle panel is the same simulation, but at four orbits, and the right panel is from PEnGUIn at four orbits. At one orbit, FARGO3D finds a structured flow pattern inside the planet's Hill sphere, with supersonic flow directed toward the planet from above. This pattern breaks down to a flow resembling turbulence at four orbits. In contrast, PEnGUIn sustains a structured flow pattern for the entire duration of the simulation.

<figure><img src="image_11.png" /></figure>

**Figure 9.** Vertical velocity $v_z$ plots, same as Figure 8, but at an azimuthal slice offset from the planet's position by 0.01 radian in azimuth. Also note that the x-axis is shifted inward from the planet's orbit. We plot the $\{r_s/r_H, \gamma\} = \{0.03, 1.2\}$ model at high resolution, obtained at four orbits for both FARGO3D (left) and PEnGUIn (right). We find meridional vortices, which are illustrated here by the strong, localized velocity shear, traveling downstream from the planet in both codes.

the torque in these last 5 orbits. In the left panel, we see that the final torque in high resolution converges to a moderately larger value than in low resolution. In the right panel, we plot the torque from the high-resolution run with two different values of the smoothing length $r_s$. As was the case with our isentropic results, there is no clear dependence on the smoothing length. Inserting our disk model into the formula of Masset & Casoli (2010) and using an updated estimate of the width of the horseshoe region (Lega et al. 2015; Masset & Benítez-Llambay 2016) gives a torque value of 0.34 $\Gamma_0$, in rough agreement with our results. We note that, when considering radiative effects, the corotation torque can be positive and possibly dominate over the negative Lindblad torque (Masset & Casoli 2009, 2010; Paardekooper et al. 2010, 2011), leading to outward migration such as we observe here.

## 3.2. Flow Pattern Around the Planet

In this section, we investigate the connection between the flow pattern around the planet and the planetary torque. Previous simulations of embedded planets have seen vertical flow toward

the planet's poles (Kley et al. 2001; Klahr & Kley 2006; Tanigawa et al. 2012; Morbidelli et al. 2014; Szulágyi et al. 2014; Ormel et al. 2015; and FAW15), and FAW15 further showed that this results in the generation of meridional vortices. A similar flow pattern is recovered in some, but not all, of our simulations. In particular, we find significant disagreement between FAR-GO3D and PEnGUIn simulations.

Figure 8 plots the vertical velocity $v_z$ in a meridional slice across the planet's location, corresponding to the $\{r_s/r_H, \gamma\} = \{0.03, 1.2\}$, high-resolution model. It shows that, despite the close agreement in torque measurements between FARGO3D and PEnGUIn, their flow patterns are different. At four orbits, PEnGUIn shows an organized flow structure where a fast ($\sim 2 c_0$) vertical downflow is directly above the planet, and two vortices are generated at $z \sim 0.4 r_H$, qualitatively similar to those reported by FAW15 (e.g., their Figure 11). FARGO3D shows a similar pattern only at $t \sim 1$ orbit; for $t \ge 2$ orbits, it displays a turbulent flow within the planet's Hill sphere, with a maximum speed of $\sim 0.2 c_0$. This disagreement between the two codes is generally present when $r_s = 0.03 r_H$.
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Vertical velocity $v_z$ plots, same as Figure 8, but with different disk models. The left panel is identical to the right panel of Figure 8, shown here for comparison. The middle panel plots a model with a larger $r_s$: {$r_s/r_H, \gamma$} = {0.05, 1.2}, and the right panel with a larger $\gamma$: {$r_s/r_H, \gamma$} = {0.03, 1.4}. All three panels are from PEnGUIn simulations. The maximum vertical flow speed on each panels are 2.0, 0.16, and 0.28 $c_0$ from left to right, respectively, indicating that flow speeds are slower when $r_s$ or $\gamma$ is larger.</figcaption></figure>

and at high resolution. It is unclear to us which of these flow patterns, turbulent or organized, is correct, but we are encouraged by the fact that the torque measurements appear to be robust against these discrepancies.

Despite their differences, both FARGO3D and PEnGUIn find that the planet sheds meridional vortices downstream, in agreement with FAW15. Figure 9 again plots $v_z$, but in a meridional slice offset from the planet's location by 0.01 radian. The vortex, seen as a strong, localized shear in $v_z$, is evident in both codes, with the vortex in FARGO3D being noticeably weaker.

Figure 10 uses three models from PEnGUIn to show how the flow pattern depends on $r_s$ and $\gamma$. The right panel is identical to the left panel of Figure 8; the middle panel shows the flow when $r_s$ is larger; and the left panel when $\gamma$ is larger. In general, the flow pattern is preserved when varying $r_s$ and $\gamma$, and only the flow speed changes. Increasing $r_s$ or $\gamma$ both reduces the flow speed; the same applies for FARGO3D simulations as well.

## 4. Conclusion and Discussion

We perform hydrodynamical simulations of disk-planet interaction using two different codes, FARGO3D and PEnGUIn, and measure the torque on a planet embedded in a 3D isentropic disk. We find that the torque is independent of the smoothing length $r_s$ (Figure 2), and only weakly dependent on the adiabatic index $\gamma$ (Figure 3). In order to obtain convergence on these measurements, we vary the resolution of our simulations, and find that convergence requires at least three cells per smoothing length $r_s$ when $\gamma \gtrsim 1.4$, and 4 ~ 5 cells when $\gamma = 1.2$. Overall, our two codes show close agreement. For 5 $M_\oplus$ planets, we observe a nonlinear behavior in the torque, where a weak corotation torque is present despite our constant entropy and vortensity profiles. This corotation torque is more negative when $\gamma$ is smaller, with a magnitude as large as ~0.3 $\Gamma_0$ when $\gamma = 1.2$. For 1 $M_\oplus$ planets, the net torque is more negative with a larger $\gamma$, and no corotation torque is found. The net torque agrees with the results of Tanaka et al. (2002) when extrapolated to $\gamma = 1$.

To further assess the issue of numerical convergence, we perform additional simulations with explicit radiative transfer using FARGOCA, and again find good agreement between the measured torques and the analytic formula of Masset & Casoli

(2010). These results confirm what was previously found in lower-resolution simulations by Lega et al. (2015).

The flow field around the planet reveals some dependency on $r_s$. In PEnGUIn simulations, the flow pattern around the planet is largely preserved when $r_s$ is varied, but the flow speed decreases with a larger $r_s$. With FARGO3D, we find turbulent flow when $r_s$ is small, and calmer, slower flow when $r_s$ is larger. Increasing $\gamma$ has an effect similar to that of increasing $r_s$; it also reduces the flow speed.

Because our torque measurements are converged, both between the two codes and in resolution, we are confident in their accuracy. To account for a realistic migration rate, however, a few caveats need to be addressed. First, we measure the planetary torque after the planet has been introduced for five orbits, which is a short time compared to the migration timescale of the planet. At low resolution, we extend one of our simulations to 100 orbits, and find that the planet begins to open a gap, a consequence of our disk being inviscid. This results in a torque that strongly fluctuates in time. Numerical convergence in this case is difficult to achieve, since long-term high-resolution runs are prohibitively expensive. Second, our planet has a fixed, circular orbit. If the planet is allowed to migrate in accordance with the torque it receives, it can potentially modify the torque value. In our current setup, resolution is concentrated at the planet's initial position, so it is not fit to simulate a migrating planet. Future work, with a different numerical treatment, is required to study this aspect of migration. Finally, our parameter space is restricted to disks with an initial zero vortensity gradient. Allowing for a finite gradient is likely to produce a larger horseshoe drag that can alter the magnitude, or even the direction of migration. Although this type of study has been done previously in 2D (e.g., Casoli & Masset 2009; Masset & Casoli 2009), detailed 3D studies have only been done for highly sub-thermal planets ($M_p \ll M_{\text{thermal}}$) (Masset & Benítez-Llambay 2016). The meridional vortices observed here and in previous work may play a role in the 3D horseshoe drag for near to super-thermal planets ($M_p \gtrsim M_{\text{thermal}}$). Indeed, we have uncovered a 3D corotation torque for nearly-thermal planets in this study. Its magnitude appears small compared to the net torque, but it is possible that our particular choice of a zero vortensity gradient minimized it.

Although our torque measurements appear robust, the flow structure around the planet shows discrepancies between our
two codes. Moreover, the flow is sensitive to $r_s$, which controls the flow speeds within the Hill sphere, or even generate turbulence in FARGO3D simulations. The flow speed is important for determining how fast the gas within the Hill sphere circulates with the rest of the disk, and has implications for the thermal cooling of the planet's atmosphere (Ormel et al. 2015). Additionally, our result implies that $r_s$ can potentially affect numerically measured accretion rates, such as those by Machida et al. (2010) and Szulágyi et al. (2014). It remains to be seen whether a sufficiently small $r_s$ will produce a converged flow structure. Further code testing is also needed to understand the difference between our two codes.

We thank an anonymous referee for helpful feedback that improved this manuscript. J.F. gratefully acknowledges support from the Natural Sciences and Engineering Research Council of Canada, the Center for Integrative Planetary Science at the University of California, Berkeley, and the Sagan Fellowship Program under contract with the Jet Propulsion Laboratory (JPL) funded by NASA and executed by the NASA Exoplanet Science Institution. F.M. and D.V. gratefully acknowledge support from UNAM's DGAPA grant PAPIIT IN101616, as well as support from CONACYT grant 178377. E.L. is thankful to ANR for supporting the MOJO project (ANR-13-BS05-0003-01). This work was performed using HPC resources from GENCI [IDRIS] (Grant 2016, [i2016047233]).

## References

* Baruteau, C., & Masset, F. 2013, in Lecture Notes in Physics, Vol. 861, ed. J. Souchay, S. Mathis, & T. Tokieda (Berlin: Springer), 201

* Bate, M. R., Lubow, S. H., Ogilvie, G. I., & Miller, K. A. 2003, MNRAS, 341, 213

* Benítez-Llambay, P., Masset, F., Koenigsberger, G., & Szulágyi, J. 2015, Natur, 520, 63

* Benítez-Llambay, P., & Masset, F. S. 2016, ApJS, 223, 11

* Bitsch, B., Boley, A., & Kley, W. 2013, A&A, 550, A52

* Casoli, J., & Masset, F. S. 2009, ApJ, 703, 845

* D'Angelo, G., Kley, W., & Henning, T. 2003, ApJ, 586, 540

* D'Angelo, G., & Lubow, S. H. 2010, ApJ, 724, 730

* de Val-Borro, M., Edgar, R. G., Artymowicz, P., et al. 2006, MNRAS, 370, 529

* Dressing, C. D., & Charbonneau, D. 2015, ApJ, 807, 45

* Fung, J. 2015, PhD thesis, Univ. Toronto

* Fung, J., Artymowicz, P., & Wu, Y. 2015, ApJ, 811, 101

* Klahr, H., & Kley, W. 2006, A&A, 445, 747

* Kley, W. 1989, A&A, 208, 98

* Kley, W. 1998, A&A, 338, L37

* Kley, W., Bitsch, B., & Klahr, H. 2009, A&A, 506, 971

* Kley, W., D'Angelo, G., & Henning, T. 2001, ApJ, 547, 457

* Korycansky, D. G., & Papaloizou, J. C. B. 1996, ApJS, 105, 181

* Lega, E., Crida, A., Bitsch, B., & Morbidelli, A. 2014, MNRAS, 440, 683

* Lega, E., Morbidelli, A., Bitsch, B., Crida, A., & Szulágyi, J. 2015, MNRAS, 452, 1717

* Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321

* Machida, M. N., Kokubo, E., Inutsuka, S.-I., & Matsumoto, T. 2010, MNRAS, 405, 1227

* Masset, F. 2000, A&AS, 141, 165

* Masset, F. S., & Benítez-Llambay, P. 2016, ApJ, 817, 19

* Masset, F. S., & Casoli, J. 2009, ApJ, 703, 857

* Masset, F. S., & Casoli, J. 2010, ApJ, 723, 1393

* Masset, F. S., D'Angelo, G., & Kley, W. 2006, ApJ, 652, 730

* Mihalas, D., & Weibel Mihalas, B. 1984, Foundations of Radiation Hydrodynamics (New York: Oxford Univ. Press)

* Morbidelli, A., Szulágyi, J., Crida, A., et al. 2014, Icar, 232, 266

* Ormel, C. W., Shi, J.-M., & Kuiper, R. 2015, MNRAS, 447, 3512

* Paardekooper, S.-J., Baruteau, C., Crida, A., & Kley, W. 2010, MNRAS, 401, 1950

* Paardekooper, S.-J., Baruteau, C., & Kley, W. 2011, MNRAS, 410, 293

* Petigura, E. A., Marcy, G. W., & Howard, A. W. 2013, ApJ, 770, 69

* Silburt, A., Gaidos, E., & Wu, Y. 2015, ApJ, 799, 180

* Szulágyi, J., Morbidelli, A., Crida, A., & Masset, F. 2014, ApJ, 782, 65

* Tanaka, H., Takeuchi, T., & Ward, W. R. 2002, ApJ, 565, 1257

* Tanigawa, T., Ohtsuki, K., & Machida, M. N. 2012, ApJ, 747, 47

* Uribe, A. L., Klahr, H., Flock, M., & Henning, T. 2011, ApJ, 736, 85

* Ward, W. R. 1991, LPSC, 22, 1463