Outward Migration of Super-Jupiters

A planet embedded in a disk excites waves (Goldreich & Tremaine 1979). These waves carry angular momentum and therefore torque the disk where they damp, thereby altering the disk's surface density profile. Furthermore, the waves' angular momentum comes at the expense of the planet's, and so the planet must migrate. Goldreich & Tremaine (1980), Artymowicz (1993), and Ward (1997) worked out how much angular momentum is excited in waves, and their theory adequately predicts what is found in simulations, even for very deep gaps and Jupiter-mass planets—provided the surface density profile is known (see, e.g., the comparison in Figure 14 of DLL).

There is, however, an important missing component: how far do waves travel before damping? The answer is crucial, because even if the waves' angular momentum is known, one can only determine how the waves affect the disk's (azimuthally averaged) surface density profile if one knows where that angular momentum is deposited. The surface density profile, in turn, determines the amplitude of the excited waves, and hence the angular momentum carried by them. Ward (1997) answered this question by assuming that waves damp immediately upon being excited. With this assumption, gaps are exponentially deep for planet masses slightly greater than the gap-opening threshold (Syer & Clarke 1995; Ward 1997). Therefore, disk material is strongly prevented from crossing the planet's orbit, and the planet is locked in the disk—aka the type II paradigm.

The assumption of zero damping length is unfounded, as recognized by Ward (1997, Section IV-b) and many others (e.g., Lunine & Stevenson 1982; Goodman & Rafikov 2001; Rafikov 2002; Duffell 2015; Kanagawa et al. 2015). If the waves travel before damping, the predicted gap becomes much shallower than the exponential depth at zero damping length (Crida et al. 2006; Duffell 2015; Kanagawa et al. 2015; Ginzburg & Sari 2018; DLL). This is true even if the damping length is modest—of order a scale height. The fact that simulated gaps are not exponentially deep, and the reason for the failure of type II, is therefore simply that waves travel before damping. ⁴ In order to accurately determine where in the disk the waves damp, and thereby obtain the gap depth and the torques on the planets, one may build a theory, as has been done for low-mass planets (e.g., Goodman & Rafikov 2001). But for near-Jupiter-mass planets, which excite strongly nonlinear waves, the theory has thus far proven intractable. Therefore, we turn to hydrodynamical simulations.

The main difference between the simulations in DLL and those done by other groups is the boundary conditions. Conventionally, the disk is fixed to its initial profile at the boundaries. But that does not allow the disk to relax to a steady-state profile that is independent of the boundary locations. Instead, at our outer boundary, the disk is fed with a constant $\dot{M}$ by using the analytic solution for a steady-state disk beyond the zone where the waves have damped. That solution permits a pileup or deficit of gas beyond the planet's orbit. And at the inner boundary, material is drained without injecting angular momentum, which allows the disk there to reach the surface density profile of a free steady disk, which has $\dot{M}=3\pi \nu {\rm{\Sigma }}$. In practice, this is done by setting νΣ to be constant across the inner boundary (see also the discussion of the inner boundary condition in Dempsey et al. 2020b). Our boundary conditions are valid whether or not the type II hypothesis is correct. If it is correct, i.e., if no material crosses the planet's orbit, then the inner disk in our simulations would drain away, and the outer pileup would continually grow. Instead, we find that the disk reaches a viscous steady state, with a constant flow past the planet.

Our boundary conditions here differ from DLL in one significant way: we enforce the outer boundary condition in the barycentric frame (see Appendix B for details). If we instead use the stellocentric frame, we find that the indirect terms generate an artificial eccentricity near the outer boundary—an effect that was less pronounced in DLL because of the smaller-mass planets considered there. In many of these new simulations, the disks are found to be eccentric. But in all cases, the eccentricity near the outer boundary is sufficiently small that it is appropriate to treat orbits there as circular around the barycenter (e.g., Figure 3 below).

For most of our simulations, we place the computational boundaries at r_in = 0.3r_p and r_out = 12r_p , where r_p is the planet's orbital radius. To prevent wave reflections near the inner boundary, we place a wave-killing region between r = [0.3r_p , 0.4r_p ]. ⁵ Our method of wave killing preserves angular momentum (DLL) and so captures all of the torque injected by the planet in the computational domain. Wave killing is unnecessary at the distant outer boundary, because waves dissipate before reaching it.

We run simulations with the q, α, and h values shown in Figure 1. Each simulation is run until the time- and azimuthally averaged $\dot{M}(r)$ profile is spatially constant to within <10% (see Appendix C). In many of the simulations, ΔT does not reach a constant but instead varies quasiperiodically in a steady state (see Figure 5). As discussed below, the quasiperiodicity is associated with the disk being eccentric and apsidally precessing.

Zoom In Zoom Out Reset image size

Convergence typically takes up to seven viscous times, evaluated at r_p . The simulations are computationally expensive both because we run for a number of viscous times and because some disks become eccentric, which forces the time step to be smaller. For example, on a single P100 GPU, one viscous time at r_p takes ∼20 days of wall time for our lowest α simulations.

In Figure 1, we plot a plus or minus sign to indicate the direction of the planet's migration. At each h, the transition between inward and outward migration is roughly determined by the value of K ≡ q²/(α h⁵), with the critical value of K depending on h. We note that K is equal to the ratio of the one-sided torque in the absence of a gap (∝ q²/h³) to the viscous torque (∝ α h²; Lin & Papaloizou 1986; Duffell & MacFadyen 2013; Kanagawa et al. 2015) and found in simulations to set the gap depth when q ≲ 10⁻³ (Duffell & MacFadyen 2013; Fung et al. 2014; Kanagawa et al. 2017).

The migration rate is obtained by angular momentum conservation, $\displaystyle \frac{1}{2}{M}_{p}{{\ell }}_{p}{\dot{r}}_{p}/{r}_{p}=-{\rm{\Delta }}T$, where M_p is planet mass. We may rewrite this as

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{r}}_{{\rm{p}}}}{{r}_{{\rm{p}}}}=-\displaystyle \frac{1}{{\tau }_{\nu }^{* }}\left(\displaystyle \frac{{\rm{\Delta }}T}{\dot{M}{{\ell }}_{p}}\right),\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{1}{{\tau }_{\nu }^{* }}\equiv \displaystyle \frac{{M}_{d}}{{M}_{p}}\displaystyle \frac{1}{{\tau }_{\nu }}\end{eqnarray} \tag{ 3 }$$

is the "mass-reduced viscous rate," which we define in terms of the following quantities: the viscous time ${\tau }_{\nu }=\tfrac{2}{3}{r}_{p}^{2}/\nu ({r}_{p})$ and the proxy disk mass ${M}_{d}=4\pi {r}_{p}^{2}{{\rm{\Sigma }}}_{Z}({r}_{p})$, where ${{\rm{\Sigma }}}_{Z}=\dot{M}/(3\pi \nu )$ is the surface density profile in the absence of the planet (DLL). For comparison, the type II migration rate for planets less massive than the disk is −1/τ_ν (Ward 1997). And for planets more massive than the disk (upon which we focus), type II predicts the rate to be $-{({M}_{d}/{M}_{p})}^{\gamma }/{\tau }_{\nu }^{* }$, where γ is a positive number whose value depends on the assumed background disk profile (Syer & Clarke 1995; Ivanov et al. 1999).

Figure 2 (top left panel) shows the migration rates from the simulations. The curves in the top left panel are from DLL, in which we ran simulations with q ≤ 10⁻³ and h = 0.05, and fit the resulting ΔT to a power-law expression. The extrapolation of those curves to higher q is clearly discrepant with our simulations here.

Zoom In Zoom Out Reset image size

For our present simulations with higher-mass planets, we see that at high enough q, planets migrate outward. For example, focusing first on the simulation set with h = 0.05 (filled diamonds), we see that the subset with α = 0.01 transitions to outward migration at q ≳ 4 × 10⁻³. At lower α (and still h = 0.05), the transition to outward migration occurs at a lower q.

The transition from inward to outward migration happens alongside another transition: the disk becomes eccentric. In the bottom left panel of Figure 2, we plot the maximum eccentricity in the disk ${e}_{\max }$ versus q, where the eccentricity is measured by averaging the eccentricity vector within each ring of the disk and then time-averaging the resulting (scalar) eccentricity (Teyssandier & Ogilvie 2017). Comparing each subset of simulations in the top left and bottom left panels, we see that migration is inward in nearly circular disks and outward in eccentric ones (i.e., when ${e}_{\max }\gtrsim 0.1$).

As seen in Figure 2, both transitions occur above a critical q that depends on α and h. We find empirically that the transition occurs at roughly a fixed value of

$$\begin{eqnarray}&&K^{\prime} \equiv {{Kh}}^{2}=\displaystyle \frac{{q}^{2}}{\alpha {h}^{3}}.\end{eqnarray} \tag{ 4 }$$

The right panels of Figure 2 replot the data in the left panels, but now with $K^{\prime}$ on the x-axis. We observe that at $K^{\prime} \gtrsim 20$, i.e., at

$$\begin{eqnarray}&&q\gtrsim 1.5\times {10}^{-3}{\left(\displaystyle \frac{\alpha }{0.001}\right)}^{1/2}{\left(\displaystyle \frac{h}{0.05}\right)}^{3/2},\end{eqnarray} \tag{ 5 }$$

the migration rate transitions from inward to outward, and the disk eccentricity transitions from ≲0.1 to ≳0.2. The parameter $K^{\prime}$ has been found empirically to correlate with the width of the gap for planets not massive enough to excite a significant disk eccentricity (Kanagawa et al. 2016; DLL). Kley & Dirksen (2006) and Teyssandier & Ogilvie (2017) found that for α = 0.004 and h = 0.05, disks become eccentric when q ≳ 3 × 10⁻³, in agreement with Equation (5).

The theory laid out in Lubow et al. (1999) and Teyssandier & Ogilvie (2016, 2017) shows that the disk's eccentricity excitation is sensitive to the density profile, because the density profile controls whether the resonances that excite eccentricity are stronger than those that damp it. Hence, the fact that $K^{\prime}$—and hence the gap width—appears to control whether disks are eccentric or not is not too surprising. Nonetheless, although the theory of Teyssandier & Ogilvie (2016) successfully predicts many aspects of their simulations, it does not reproduce the correct threshold for where eccentricity is excited, so they may be missing an effect.

Figure 3 shows the steady-state profiles of surface density and eccentricity for simulations that have h = 0.05 and $K^{\prime}$ near the migration transition. In the top row, migration is inward and the disk is nearly circular; in the bottom, migration is outward and the disk is eccentric. The profiles are time-averaged in steady state over 10,000 orbits, which is much longer than the precessional times. As seen in the Σ profiles, higher values of $K^{\prime}$ systematically result in deeper gaps.

Zoom In Zoom Out Reset image size

In the circular disks, the gaps are nearly symmetric relative to the planet's position, while in the eccentric disks, the outer half of the gap becomes significantly wider than the inner one. In addition, the circular disks have a modest density pileup outside of the planet's orbit, and the eccentric ones have a deficit. The sign and magnitude of the pileup is dictated by the value of ΔT (as shown in Equation (19) of DLL, which follows from Equation (7) below).

Turning to the eccentricity profiles, we see that above the transitional $K^{\prime}$, it is the outer disk that develops a significant eccentricity. And, although not shown in the figure, the eccentric disks precess coherently, ⁶ and the amplitude of the eccentricity remains unchanged over multiple viscous times. This points to the eccentricity being a free mode of the disk (Teyssandier & Ogilvie 2016; Lee et al. 2019), with an excitation that balances viscous damping to maintain the mode in a steady state (as in circumbinary disks; see Muñoz & Lithwick 2020).

The right panels show T_dep, the cumulative torque deposited into the disk (within each r) by the damping of waves. To calculate T_dep, we evaluate the angular momentum flux carried by waves and then remove from that the gravitational torque excited by the planet (DLL). We see that throughout most of the outer disk, T_dep increases with r, showing that the torque deposition per unit r is positive there; i.e., that region is responsible for pushing the planet inward. Similarly, most of the inner disk pushes the planet outward, as T_dep decreases with r when r < r_p . Far beyond the planet—outside a few × r_p —T_dep flattens and approaches ΔT; i.e., all of the torque excited by the planet is eventually deposited in the disk. We henceforth use the T_dep profiles to examine what is causing ΔT to become negative in eccentric disks.

Figure 3 shows that the torque on the inner disk is always

$$\begin{eqnarray}&&{T}_{\mathrm{dep}}(r={r}_{p})\approx -\dot{M}{{\ell }}_{p}.\end{eqnarray} \tag{ 6 }$$

Equation (6) is a generic result for deep gaps in steady-state disks. It follows from the conservation of angular momentum flux, which we approximate here as

$$\begin{eqnarray}&&{T}_{\mathrm{dep}}(r)\approx 3\pi \nu \langle {\rm{\Sigma }}\rangle {\ell }-\dot{M}{\ell }\end{eqnarray} \tag{ 7 }$$

(see DLL for the full expression). Therefore, if the gap is sufficiently deep, the first term on the right-hand side may be neglected at r_p , confirming Equation (6). In the right panels, we also convert the Σ profiles to T_dep via Equation (7). The result is shown as dotted curves. The fact that these agree with the directly calculated T_dep confirms that the disk is in viscous steady state.

Given the value of the torque on the inner disk (Equation (6)), it must be that ΔT becomes negative in eccentric disks because the torque on the outer disk weakens. But why? From the theory in Section 1.1, we must examine torque excitation and deposition.

• 1.  
  Excitation. For a given Σ profile, eccentricity could lower ΔT if it lowers the amplitude of the excited waves in the outer disk and hence lowers the angular momentum carried by them. (We also include in this category the change in co-orbital torques; see below.)
• 2.  
  Deposition. For a given ΔT, eccentricity could increase the distance waves travel before they damp. A longer damping length means a wider gap (from Equation (7)), which implies that at the location where the waves are excited, Σ is lower, and hence the wave amplitude is lower.

We examine item 1 in Figure 4 by comparing the simulated values of ΔT with a linear calculation of ΔT in a circular disk whose Σ(r) profile is the same as that from the simulation ⁷ (see DLL for our linear calculation method.) We see that for the circular disks ($K^{\prime} \lesssim 20$), the linear calculation adequately predicts what is found in the simulations. But for the eccentric disks, it overpredicts ΔT. We infer that the disk being eccentric lowers ΔT, as suspected in item 1. A more detailed examination shows that much of the lowering is indeed due to outer waves carrying less angular momentum in the simulated eccentric disk (relative to the linear calculation). But there is an additional effect: the co-orbital torques become stronger (more negative) in the eccentric disks. We find this effect to be subdominant up to $K^{\prime} \sim 50;$ for $K^{\prime}$ larger than that, the two effects are comparable.

Zoom In Zoom Out Reset image size

Turning to item 2, the right panels of Figure 3 show that the damping length is indeed longer in the eccentric disks. Unfortunately, we have not been able to quantify the relative importance of items 1 and 2. For a partial quantification of the latter, the curves in Figure 4 show the extrapolations from the circular planet simulations (as in Figure 2). These lie above the open diamonds for the eccentric disks ($K^{\prime} \gtrsim 10$), and one of the reasons for the discrepancy is that the gap shape for the open diamonds is wider than would be implied in the extrapolations. In truth, explanations 1 and 2 are entangled; e.g., when waves are excited to lower amplitudes, they travel further before becoming nonlinear and damping. But for a first-principles derivation, one must extend the theory of Goodman & Rafikov (2001) to investigate more massive planets and viscous disks.

Figure 2 encapsulates our main result: planet migration transitions from inward to outward when the planet exceeds around 2 M_J for the fiducial values α = 0.001 and h = 0.05 (Equation (5)). This transition coincides with the disk becoming eccentric. In Section 3.2, we showed why eccentric disks lead to outward migration.

We find that the total torques are always close in magnitude to the advection torque ($| {\rm{\Delta }}T| \sim \dot{M}{{\ell }}_{p}$), whether migration is inward or outward (Figure 2). As a result, the migration rates are comparable in magnitude to the "mass-reduced viscous rate" (Equation (2)). ⁸ It is not surprising that $| {\rm{\Delta }}T| \sim \dot{M}{{\ell }}_{p}$ when migration is outward, because the inner torque is always equal to $-\dot{M}{{\ell }}_{p}$ for deep gaps (Section 3.2); so, if the outer torque is subdominant, the total torque will also be close to $-\dot{M}{{\ell }}_{p}$. But it is surprising that, even when migration is inward, it is never the case that ${\rm{\Delta }}T\gg \dot{M}{{\ell }}_{p}$. We speculate that this is because a planet always acts as a leaky sieve, so the pileup beyond the planet can never be too big. That would, in turn, imply that the normalized torque cannot be too large because of the close connection between the pileup and the normalized torque (DLL). But why the planet should always act as a leaky sieve remains an open question.

• 1.  
  α-model. Perhaps our most questionable assumption is that of α-viscosity. If angular momentum is transported by a nonviscous mechanism, such as disk winds, then the gap shapes would be different, and hence so would the torques.
• 2.  
  Low-mass disk. Fixing the planet's orbit is a good approximation when the planet's migration timescale ($\sim {\tau }_{\nu }^{* }$) is shorter than the gas radial drift timescale (τ_ν ) or, equivalently, when M_d ≲ M_p .
• 3.  
  Circular planet orbit. An eccentric planet will introduce a forced component to the disk eccentricity on top of the free component that is excited above the transitional mass. The planet's eccentricity will also modify the torque on the disk. D'Angelo et al. (2006) included the backreaction of the disk on the planet and found that the planet's eccentricity is rapidly excited, which affects the migration direction. But their inner disks are severely depleted, perhaps due to their inner boundary condition, an insufficient run time compared to the viscous time, or material being lost onto the planet. Moreover, their disks are at least as massive as the planet, invalidating our low-mass disk assumption. In the future, we plan to determine whether the planet's eccentricity is excited or damped in low-mass steady disks, as has been done for stellar binaries (e.g., Muñoz et al. 2019; D'Orazio & Duffell 2021). In the latter case, it has been found that the binary's eccentricity damps if it is below a threshold eccentricity of around 0.1, even though the disk is very eccentric.
• 4.  
  2D disk. For high-mass planets, gaps are much wider than the disk scale height, so the 2D treatment should adequately capture the excitation and deposition torques.
• 5.  
  No accretion onto planet. Our large softening parameter prevented planetary accretion. Realistic modeling of accretion requires high resolution, in addition to a more accurate treatment of the thermal state of the gas, disk self-gravity, and three dimensions (e.g., Fung et al. 2019). If a significant fraction of the disk's radial mass flow ends up in the planet, the migration rates would change significantly. Additionally, the torque from the circumplanetary region is unrealistic. But we have checked that this torque is small, typically contributing between 0.1 and 0.25 to the dimensionless migration rates in Figure 2.
• 6.  
  Locally isothermal equation of state. Adopting a finite cooling time, rather than the instantaneous cooling that we assume, would affect the deposition profile. That is because linear waves deposit angular momentum into the disk even in the absence of viscous dissipation, and the amount depends on the cooling time (Miranda & Rafikov 2020). We suspect that this dependence is minimal for the very massive planets considered in this paper, for which viscous dissipation (including by shocks) likely dominates over the aforementioned linear wave deposition.
• 7.  
  Inner boundary condition. For our results not to depend on the inner boundary, the inner wave-killing zone must be far enough from the planet that most of the torque excitation occurs within the simulation domain (DLL). We have checked that by examining the excitation profile. In addition, we have tested roughly half of our h < 0.1 simulations by continuing them in a domain with a smaller wave-killing boundary (down to 0.2r_p ) and found in these cases that the change in torque was small. Note that the h = 0.1 simulations require a larger radial domain (see footnote ²) because the wave excitation profile is broader at larger h.

Previous studies have found outward migration of massive planets, but these require either steeply falling surface density profiles (Chen et al. 2020) or extremely massive disks to induce either type III migration (Masset & Papaloizou 2003) or gravitational instability (Lin & Papaloizou 2012; Cloutier & Lin 2013). By contrast, we find outward migration more generally for super-Jupiters, subject to the assumptions listed above.

Recently, Duffell et al. (2020) examined a setup similar to ours and found that ΔT reverses sign in the brown dwarf regime, but for h = 0.1 and α ≥ 0.03. Using Equation (5) for their parameters, we find general agreement with their transitional masses (see their Figure 5), suggesting that the torque reversal reported by Duffell et al. (2020) is equivalent to ours.

Zoom In Zoom Out Reset image size

The 2 M_J threshold that we have uncovered is intriguingly similar to the lowest-mass giant planets discovered via direct-imaging techniques (e.g., Bowler & Nielsen 2018). Due to their large masses and wide separations, these planets were once thought to originate from gravitational instabilities. But recent evidence points to a mass distribution that may be consistent with the core accretion scenario (Nielsen et al. 2019; Wagner et al. 2019). Thus, those super-Jupiters detected via direct imaging could be the outcome of core accretion but followed by outward migration in an eccentric disk. This scenario would not require the presence of additional giants, as would be the case for jointly migrating planets trapped in mean-motion resonance (Crida et al. 2009). It has indeed been speculated that outward migration could explain the properties of directly imaged planets, e.g., as pointed out recently by Bohn et al. (2021) in the context of their discovery of an ∼6 M_J planet at ≳115 au from its star. We propose that such planets could have migrated via interaction with an eccentric disk.