Primordial Giant Planet Obliquity Driven by a Circumplanetary Disk

The rotation state of planets is one of the basic observables that offers information about planet formation and planetary system evolution. For giant planets, which open a gap in the protoplanetary disk (Lin & Papaloizou 1986), the primordial spin state is determined by the interaction between the growing planet and its circumplanetary disk, which forms because the planet is much smaller than the Hill radius (Artymowicz & Lubow 1996; Lubow et al. 1999; D'Angelo et al. 2002). The magnitude of the planetary angular momentum that results from circumplanetary disk accretion depends upon whether the planet–disk interaction is hydrodynamic (Dong et al. 2020) or magnetospheric (Koenigl 1991; Batygin 2018), but once set is expected to be conserved in the absence of tidal effects. In the solar system Jupiter and Saturn rotate rapidly, but at a rate that is still only 30%–40% of their nominal breakup angular velocity. For a small sample of very massive extrasolar planets, Bryan et al. (2018) inferred somewhat lower but nonzero rotation rates.

Misalignment between the spin and orbital angular momentum vectors of planets—the planetary obliquity—is common in the solar system. Jupiter's obliquity is low (3°), but Saturn (27°) and the ice giants have highly significant obliquities. Measurement of extrasolar planetary obliquities is challenging, but there are some indications of nonzero obliquity for the directly imaged planetary mass companion in the 2M0122 system (Bryan et al. 2020). Unlike the magnitude of the angular momentum, the obliquity is relatively fragile to late-time evolutionary changes. Resonances between the precession frequencies of planetary spin axes and orbits affect both terrestrial and giant planets in the solar system (Laskar & Robutel 1993; Ward & Hamilton 2004; Ward & Canup 2006). A related process, where the orbital precession is driven by the gravitational potential of the protoplanetary disks, can generate early-time obliquity (Millholland & Batygin 2019). Giant impacts (Safronov 1966; Benz et al. 1989; Morbidelli et al. 2012) and planet–planet scattering events (Li 2021; Hong et al. 2021) can also lead to planet spin–orbit misalignment.

In this Letter, we investigate another mechanism that can generate primordial obliquity, as a byproduct of planetary interaction with a tilted circumplanetary disk. The physical properties of circumplanetary disks, specifically their large size relative to the Hill sphere, and moderately large geometric thickness, favor growth of the tilt instability of disks in binary systems (Lubow 1992; Lubow & Ogilvie 2000; Martin et al. 2020). The tilt may be able to grow large enough to excite Kodai–Lidov oscillations (Kozai 1962; Lidov 1962), at which point eccentricity is induced and accretion becomes episodic (Martin et al. 2014, 2021). Here, we show that a circumplanetary disk, whose tilt is being excited by an external tidal potential, can couple to the planetary spin via oblateness, and generate obliquity. The planet exerts a torque on the disk (e.g Tremaine & Davis 2014; Speedie & Zanazzi 2020) and an equal and opposite torque is felt by the planet. This can change the spin axis of the planet depending upon the ratio of the disk and planet-spin angular momenta. An analogous process can realign a spinning black hole with a misaligned disk (Scheuer & Feiler 1996; Martin et al. 2007; Perego et al. 2009). We show that the torque from the planet causes the planet spin and the disk to precess around the sum of their angular momenta. They align toward each other on a relatively short timescale. On a longer timescale, both can have increasing tilts as a result of the tidal tilt instability, if the disk mass is sufficiently large.

We model a circumplanetary disk around a planet of mass ${M}_{{\rm{p}}}$ that orbits a star of mass ${M}_{{\rm{s}}}$ at orbital separation ${a}_{{\rm{p}}}$. The angular velocity of the planet–star system is ${{\boldsymbol{\Omega }}}_{{\rm{p}}}={{\rm{\Omega }}}_{{\rm{p}}}{{\boldsymbol{e}}}_{z}$ where the angular frequency is ${{\rm{\Omega }}}_{{\rm{p}}}=\sqrt{G({M}_{{\rm{s}}}+{M}_{{\rm{p}}})/{a}_{{\rm{p}}}^{3}}$. The disk has surface density ${\rm{\Sigma }}(r)$ (that is fixed in time) and a unit vector to describe the tilt at each radius, ${{\boldsymbol{l}}}_{{\rm{d}}}(r,t)=({l}_{{\rm{d}},x},{l}_{{\rm{d}},y},{l}_{{\rm{d}},z})$. The tilt is assumed to be small so that ${l}_{{\rm{d}},z}\approx 1$ and ${l}_{{\rm{d}},x}$, ${l}_{{\rm{d}},y}\ll 1$. In order to explore the evolution to high inclination we would need to use hydrodynamical simulations.

A circumplanetary disk is expected to be in the wave-like disk regime (Papaloizou & Pringle 1983) since the Shakura & Sunyaev (1973) α viscosity parameter is much smaller than the disk aspect ratio H/r (e.g., Martin et al. 2019), where H is the disk scale height and r is the spherical distance from the planet. Thus, we solve the 1D wave-like warped disk equations in the frame of the planet–star system. The disk is in Keplerian rotation around the planet at all radii with angular frequency given by ${\rm{\Omega }}=\sqrt{{{GM}}_{{\rm{p}}}/{r}^{3}}$. We solve Equations (12) and (13) in Lubow & Ogilvie (2000) with an additional torque on the disk from the planet. These can be written as

$$\begin{eqnarray}&&{\rm{\Sigma }}{r}^{2}{\rm{\Omega }}\left(\displaystyle \frac{\partial {{\boldsymbol{l}}}_{{\rm{d}}}}{\partial t}+{{\boldsymbol{\Omega }}}_{{\rm{p}}}\times {{\boldsymbol{l}}}_{{\rm{d}}}\right)=\displaystyle \frac{1}{r}\displaystyle \frac{\partial {\boldsymbol{G}}}{\partial r}+{{\boldsymbol{T}}}_{{\rm{s}}}+{{\boldsymbol{T}}}_{{\rm{p}}}\end{eqnarray} \tag{ 1 }$$

and

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{G}}}{\partial t}+{{\boldsymbol{\Omega }}}_{{\rm{p}}}\times {\boldsymbol{G}}+\alpha {\rm{\Omega }}{\boldsymbol{G}}=\displaystyle \frac{{\rm{\Sigma }}{H}^{2}{r}^{3}{{\rm{\Omega }}}^{3}}{4}\displaystyle \frac{\partial {{\boldsymbol{l}}}_{{\rm{d}}}}{\partial r},\end{eqnarray} \tag{ 2 }$$

where $2\pi {\boldsymbol{G}}$ is the internal disk torque. We only need to solve the first two components of these vector equations since the disk angular momentum vector is a unit vector. We use the boundary conditions ${\boldsymbol{G}}=0$ and $\partial {{\boldsymbol{l}}}_{{\rm{d}}}/\partial r=0$ at $r={r}_{\mathrm{in}}$ and $r={r}_{\mathrm{out}}$. The initial conditions are chosen to be ${\boldsymbol{G}}=0$ and the disk is initially flat but tilted by 10° with ${l}_{{\rm{d}},x}=\sin (10^\circ )$ and ${l}_{{\rm{d}},y}=0$.

We work in a frame where the star is on the positive x-axis and so the stellar torque on the disk per unit area is ${{\boldsymbol{T}}}_{{\rm{s}}}=({T}_{{\rm{s}},x},{T}_{{\rm{s}},y})$, where ${T}_{{\rm{s}},x}=0$ and

$$\begin{eqnarray}&&{T}_{{\rm{s}},{\rm{y}}}=-\displaystyle \frac{3{{GM}}_{{\rm{s}}}}{2{a}_{{\rm{p}}}^{3}}{\rm{\Sigma }}{r}^{2}{l}_{{\rm{d}},x}.\end{eqnarray} \tag{ 3 }$$

The torque per unit area from a spinning oblate planet on the disk is given in linear theory by

$$\begin{eqnarray}&&{T}_{{\rm{p}},x}=T({l}_{{\rm{d}},y}-{l}_{{\rm{p}},y})\end{eqnarray} \tag{ 4 }$$

and

$$\begin{eqnarray}&&{T}_{{\rm{p}},y}=T(-{l}_{{\rm{d}},x}+{l}_{{\rm{p}},x}),\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&T=\displaystyle \frac{3{{GM}}_{{\rm{p}}}{r}_{{\rm{p}}}^{2}{J}_{2}}{2{r}^{3}}{\rm{\Sigma }}\end{eqnarray} \tag{ 6 }$$

(e.g., Tremaine & Davis 2014; Speedie & Zanazzi 2020). Here ${r}_{{\rm{p}}}$ is the radius of the planet and J₂ is the quadrupole gravitational harmonic. In the solar system, the giant planets Jupiter, Saturn, Uranus, and Neptune have ${J}_{2}=0.015,0.016,0.003$, and 0.004 respectively.

The planet-spin angular momentum vector is ${J}_{\mathrm{planet}}{{\boldsymbol{l}}}_{{\rm{p}}}$, where ${{\boldsymbol{l}}}_{{\rm{p}}}$ is a unit vector and

$$\begin{eqnarray}&&{J}_{\mathrm{planet}}=\displaystyle \frac{2\pi {{kM}}_{{\rm{p}}}{r}_{{\rm{p}}}^{2}}{{P}_{\mathrm{rot}}},\end{eqnarray} \tag{ 7 }$$

where k = 0.205 (appropriate if the planet is a polytrope with n = 1.5). Over the timescales considered in this work, we ignore the effects of accretion onto the planet and assume that the magnitude of the planet-spin angular momentum does not change. The direction of the spin changes as a result of the torque from the disk on the planet according to

$$\begin{eqnarray}&&{J}_{\mathrm{planet}}\left(\displaystyle \frac{d{{\boldsymbol{l}}}_{{\rm{p}}}}{{dt}}+{{\boldsymbol{\Omega }}}_{{\rm{p}}}\times {{\boldsymbol{l}}}_{{\rm{p}}}\right)=-{\int }_{{r}_{\mathrm{in}}}^{{r}_{\mathrm{out}}}2\pi r{{\boldsymbol{T}}}_{{\rm{p}}}\,{dr}.\end{eqnarray} \tag{ 8 }$$

We take the initial condition that the planet spin is aligned with its orbit with ${{\boldsymbol{l}}}_{{\rm{p}}}={\bf{0}}$.

We solve the coupled integro-differential Equations (1), (2), and (8) as an initial value problem for ${{\boldsymbol{l}}}_{{\rm{d}}}$, ${\boldsymbol{G}}$, and ${{\boldsymbol{l}}}_{{\rm{p}}}$ using finite differences. The method is first-order explicit in time. The grid extends from inner radius ${r}_{\mathrm{in}}$ up to ${r}_{\mathrm{out}}$. In order to obtain convergence, we use a nonuniform grid. There are 25 grid points spaced linearly in $\mathrm{log}(r)$ from the inner radius up to ${r}_{\mathrm{change}}=0.02\,{r}_{{\rm{H}}}$. There are 75 points that are linearly distributed with radius from ${r}_{\mathrm{change}}$ up to ${r}_{\mathrm{out}}$.

The inclinations of the disk and planet are determined via

$$\begin{eqnarray}&&i={\cos }^{-1}\left(\sqrt{1-{l}_{x}^{2}-{l}_{y}^{2}}\right)\end{eqnarray} \tag{ 9 }$$

and the phase angles as

$$\begin{eqnarray}&&\phi ={\tan }^{-1}\left(\displaystyle \frac{{l}_{y}}{{l}_{x}}\right).\end{eqnarray} \tag{ 10 }$$

For the disk, we calculate the angular momentum weighted vector components. We transform these back into the inertial frame in our plots by rotating about the z-axis by ${{\rm{\Omega }}}_{{\rm{p}}}t$.

We model a Jupiter mass and radius planet orbiting at ${a}_{{\rm{p}}}=5.2\,\mathrm{au}$. The circumplanetary disk has surface density profile ${\rm{\Sigma }}\propto {r}^{-3/2}$ distributed between the inner disk radius ${r}_{\mathrm{in}}=2\,{r}_{{\rm{p}}}$ and the outer radius ${r}_{\mathrm{out}}=0.4\,{r}_{{\rm{H}}}$, where ${r}_{{\rm{H}}}$ is the Hill radius. The inner disk radius is close to the peak in the surface density of a steady-state disk that joins onto the planet at radius ${r}_{{\rm{p}}}$. The outer radius is chosen to be close to the tidal truncation radius of a circumplanetary disk (Martin & Lubow 2011). ⁴ We take the disk parameters $\alpha =0.01$ and $H/r=0.14$. With this disk aspect ratio the disk is strongly unstable to tilting. As shown in Martin et al. (2020), the combination of the disk aspect ratio and the disk size determines the instability. While it appears that instability requires specific parameters in the analytic model, hydrodynamic simulations show that that a disk with a more realistic surface density profile is also unstable to tilting.

Figure 1 shows the evolution of the inclination and nodal phase angle for a circumplanetary disk in the limit where there is no planet torque, ${J}_{2}=0$. The blue lines show the disk angular momentum and the red lines show the planet-spin axis. The disk evolution is similar to that found in hydrodynamical simulations (Martin et al. 2020, 2021), confirming that the one-dimensional wave-like linear model is at least a fair approximation to the full disk dynamics. The planet spin does not evolve in this case as there is no torque on the planet. This result for the disk evolution is independent of the disk mass.

Zoom In Zoom Out Reset image size

We now consider two values for the quadrupole gravitational harmonic, J₂ = 0.01 and J₂ = 0.02. These J₂ values bracket those of Jupiter and Saturn, and are appropriate for planets whose rotation is rapid but well below the breakup value. The planet is assumed to rotate with period ${P}_{\mathrm{rot}}=10\,\mathrm{hr}$. Figure 2 shows the disk and planet evolution for two disk masses, ${M}_{{\rm{d}}}=0.01\,{M}_{{\rm{p}}}$ (upper panels) and ${M}_{{\rm{d}}}=0.001\,{M}_{{\rm{p}}}$ (lower panels), for both choices of J₂. The planet and the disk precess around each other with damped inclination relative to each other. As these oscillations damp out, the stellar torque becomes the dominant effect. For the high-mass disk, the planet and the disk both increase in inclination over time. However, for the low-mass disks, the planet torque dominates the stellar torque and the disk and the planet do not increase significantly in inclination over the timescales shown here. There is a superposition of two modes: a damped mode in which the disk and the planet have oppositely directed tilts and a growing mode in which they have nearly equal tilts. In each case, the latter mode has a positive growth rate and eventually dominates but the growth rate is diminished by the coupling between the planet and the disk. Even in the low-mass disk case, the planet obliquity increases in time over longer timescales than shown here.

Zoom In Zoom Out Reset image size

The timescale on which the inclination of the planet–disk system changes is longer than that of a disk without a planetary torque, but it is still relatively short, of the order of ${10}^{2}\ {P}_{\mathrm{orb}}$. We do note that a high angular momentum disk is required to change the planet spin. Significant effects are therefore more likely to occur in planets that are farther from their host star since the circumplanetary disk is larger and therefore can have more angular momentum.

In order to make some analytic estimates we now consider the timescale on which the planet-spin axis changes as a result of the disk torque only in the absence of the star. We assume that the disk is flat but tilted so that ${{\boldsymbol{l}}}_{{\rm{d}}}={{\boldsymbol{l}}}_{{\rm{d}}}(t)$. The spin angular momentum of the planet is

$$\begin{eqnarray}&&{J}_{\mathrm{planet}}=2.5\times {10}^{-7}\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{10}^{-3}\,{M}_{\odot }}\right){\left(\displaystyle \frac{{r}_{{\rm{p}}}}{{r}_{{\rm{J}}}}\right)}^{2}{\left(\displaystyle \frac{{P}_{\mathrm{rot}}}{10\mathrm{hr}}\right)}^{-1}{M}_{\odot }\,{\mathrm{au}}^{2}\,{\mathrm{yr}}^{-1},\end{eqnarray} \tag{ 11 }$$

where ${r}_{{\rm{J}}}$ is the radius of Jupiter. The angular momentum of the disk is

$$\begin{eqnarray}&&{J}_{\mathrm{disk}}={\int }_{{r}_{\mathrm{in}}}^{{r}_{\mathrm{out}}}2\pi r{\rm{\Sigma }}({r}^{2}{\rm{\Omega }}){dr}\end{eqnarray} \tag{ 12 }$$

and for typical parameters this is

$$\begin{eqnarray}\begin{array}{rcl}{J}_{\mathrm{disk}} & = & 3.8\times {10}^{-7}{\left(\displaystyle \frac{{a}_{{\rm{p}}}}{5.2\mathrm{au}}\right)}^{1/2}\left(\displaystyle \frac{{M}_{{\rm{d}}}}{{10}^{-5}\,{M}_{\odot }}\right)\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{10}^{-3}{{\rm{M}}}_{\odot }}\right)}^{2/3}{\left(\displaystyle \frac{{M}_{{\rm{s}}}}{{{\rm{M}}}_{\odot }}\right)}^{-1/6}{M}_{\odot }\,{\mathrm{au}}^{2}\,{\mathrm{yr}}^{-1}.\end{array}\end{eqnarray} \tag{ 13 }$$

Note that this calculation depends on the surface density profile. We have taken a steep profile of ${\rm{\Sigma }}\propto {R}^{-3/2}$. With ${\rm{\Sigma }}\propto {R}^{-1/2}$, the disk angular momentum is larger by a factor of 1.5 for the same total disk mass.

Integrating Equation (1) over the disk we have

$$\begin{eqnarray}&&{J}_{\mathrm{disk}}\displaystyle \frac{\partial {{\boldsymbol{l}}}_{{\rm{d}}}}{\partial t}={T}_{\mathrm{total}}\left({{\boldsymbol{l}}}_{{\rm{d}}}\times {{\boldsymbol{l}}}_{{\rm{p}}}\right),\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}&&{T}_{\mathrm{total}}={\int }_{{r}_{\mathrm{in}}}^{{r}_{\mathrm{out}}}2\pi {rT}\,{dr}.\end{eqnarray} \tag{ 15 }$$

We can also write Equation (8) as

$$\begin{eqnarray}&&{J}_{\mathrm{planet}}\displaystyle \frac{\partial {{\boldsymbol{l}}}_{{\rm{p}}}}{\partial t}=-{T}_{\mathrm{total}}\left({{\boldsymbol{l}}}_{{\rm{d}}}\times {{\boldsymbol{l}}}_{{\rm{p}}}\right).\end{eqnarray} \tag{ 16 }$$

5.1. Planet–Disk Precession Timescale

With Equations (14) and (16) we see that the sum of the planet and the disk angular momentum is conserved. Thus, in the absence of the star, the planet and the disk both precess in a retrograde direction about the sum of the disk and planet angular momentum vector. The timescale for the precession is

$$\begin{eqnarray}&&{t}_{\mathrm{prec}}=2\pi f\displaystyle \frac{{J}_{\mathrm{planet}}}{{T}_{\mathrm{total}}},\end{eqnarray} \tag{ 17 }$$

where

$$\begin{eqnarray}&&f=\displaystyle \frac{1}{\tfrac{{J}_{\mathrm{planet}}}{{J}_{\mathrm{disk}}}+1}.\end{eqnarray} \tag{ 18 }$$

The circumplanetary disk has a power-law surface density profile ${\rm{\Sigma }}\propto {r}^{-3/2}$ from ${r}_{\mathrm{in}}=2{r}_{{\rm{p}}}$ up to ${r}_{\mathrm{out}}=0.4\,{r}_{{\rm{H}}}$. For typical parameters we find

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{t}_{\mathrm{prec}}}{{P}_{\mathrm{orb}}} & \approx & 3.2{\left(\displaystyle \frac{{r}_{{\rm{p}}}}{{r}_{{\rm{J}}}}\right)}^{5/2}{\left(\displaystyle \frac{{J}_{2}}{0.01}\right)}^{-1}{\left(\displaystyle \frac{{a}_{{\rm{p}}}}{5.2\mathrm{au}}\right)}^{-1}{\left(\displaystyle \frac{{M}_{{\rm{d}}}}{{10}^{-5}{M}_{\odot }}\right)}^{-1}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{p}}}}{0.001{M}_{\odot }}\right)}^{1/6}{\left(\displaystyle \frac{{M}_{{\rm{s}}}}{{M}_{\odot }}\right)}^{1/3}{\left(\displaystyle \frac{{P}_{\mathrm{rot}}}{10\mathrm{hr}}\right)}^{-1}\left(\displaystyle \frac{f}{0.6}\right).\end{array}\end{eqnarray} \tag{ 19 }$$

This is in good agreement with the results shown in Figure 2, since we find the precession timescales for the different models to be 3.2 (model in the upper left panel), 1.6 (upper right), 6.9 (lower right), and $3.4\,{P}_{\mathrm{orb}}$ (lower right). We note that the factor f depends upon the ratio of the disk to planet angular momenta.

5.2. Planet–Disk Alignment Timescale

The alignment timescale between the disk and the planet spin is approximated by

$$\begin{eqnarray}&&{t}_{\mathrm{align}}=\displaystyle \frac{{\left(H/r\right)}^{2}{{\rm{\Omega }}}_{{\rm{d}}}}{\alpha {\omega }_{{\rm{p}}}^{2}}\end{eqnarray} \tag{ 20 }$$

(e.g., Bate et al. 2000; Lubow & Martin 2018), where ${{\rm{\Omega }}}_{{\rm{d}}}={\rm{\Omega }}({r}_{\mathrm{out}})$ and the precession rate is ${\omega }_{{\rm{p}}}=2\pi /{t}_{\mathrm{prec}}$. For typical parameters we find

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{t}_{\mathrm{align}}}{{P}_{\mathrm{orb}}} & = & 21.4{\left(\displaystyle \frac{H/r}{0.14}\right)}^{2}{\left(\displaystyle \frac{\alpha }{0.01}\right)}^{-1}{\left(\displaystyle \frac{{r}_{{\rm{p}}}}{{r}_{{\rm{J}}}}\right)}^{5}{\left(\displaystyle \frac{{J}_{2}}{0.01}\right)}^{-2}{\left(\displaystyle \frac{{a}_{{\rm{p}}}}{5.2\mathrm{au}}\right)}^{-2}\\ & & \times {\left(\displaystyle \frac{{M}_{{\rm{d}}}}{{10}^{-5}{M}_{\odot }}\right)}^{-2}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{0.001{M}_{\odot }}\right)}^{1/3}{\left(\displaystyle \frac{{M}_{{\rm{s}}}}{{M}_{\odot }}\right)}^{2/3}\\ & & \times {\left(\displaystyle \frac{{P}_{\mathrm{rot}}}{10\mathrm{hr}}\right)}^{-2}{\left(\displaystyle \frac{f}{0.6}\right)}^{2}.\end{array}\end{eqnarray} \tag{ 21 }$$

Since this is an exponential decay timescale, this is in rough agreement with the numerical simulations shown in Figure 2. The alignment timescales are 21.4 (upper left), 5.3 (upper right), 100.6 (lower left), and $25.2\,{P}_{\mathrm{orb}}$ (lower right). Note again that the factor f depends upon the ratio of the disk to planet angular momenta.

A circumplanetary disk that develops a tilt, due to the tidal tilt instability (Lubow 1992), provides a substantial reservoir of misaligned angular momentum even when the disk mass is small relative to that of the planet. In this Letter, we have shown that the disk misalignment can be communicated to the planet via the torques that result from planetary oblateness, leading to an efficient mechanism for generating obliquity. The mechanism will work (given our approximations) if, first, the circumplanetary disk properties (primarily the disk aspect ratio and size) allow for tilt growth. That tilt will change the planet's rotation state if, additionally, the disk mass is sufficiently large.

Empirical constraints on circumplanetary disk masses are currently weak, while theoretical models span a broad range (Lunine & Stevenson 1982; Canup & Ward 2002; Mosqueira & Estrada 2003; Lubow & Martin 2013; Batygin & Morbidelli 2020). Broadly speaking, static Minimum Mass Satellite Nebula models, those that include dead zones (Gammie 1996), and those where the disk is a decretion rather than an accretion disk, lead to large circumplanetary disk masses that could exceed the threshold we have estimated for Jovian conditions. Simple viscous disk models, in which $\alpha \sim {10}^{-2}$ and the planet is assembled in $\sim {10}^{6}\ \mathrm{yr}$, lead to subthreshold masses. Since the size of a circumplanetary disk is determined by the size of the Hill radius, which increases with distance from the star, all else being equal larger planetary semimajor axes lead to a greater likelihood of obliquity excitation. We note that the low accretion rate inferred from the Hα luminosity of PDS 70b (Y. Zhou, private communication, 2021), in the directly imaged PDS 70 system (Keppler et al. 2018), would imply a very low disk mass in the viscous scenario but could also be consistent with higher disk masses in the other models.

Our results are based upon describing the circumplanetary disk using one-dimensional equations, derived in the linear regime for wave-like warp propagation. The same tidal potential that is responsible for the tilt instability also drives spiral shocks in circumplanetary disks (Zhu et al. 2016), raising the question of whether a one-dimensional description is adequate. The tilt instability, however, is recovered in three-dimensional simulations using independent numerical methods (smoothed particle hydrodynamics and fixed-grid simulations), with a growth rate in approximate agreement with linear estimates (Martin et al. 2020). We therefore expect the current treatment to be adequate for moderate tilts, though simulations are indispensable for modeling very large misalignments. We have also ignored the effects of accretion onto the circumplanetary disk (Schulik et al. 2020; Tanigawa et al. 2012; Szulágyi et al. 2014). This is also an intrinsically three-dimensional effect, though we note that it could be approximately incorporated into our model using methods analogous to those used in dwarf nova disk models (Bath & Pringle 1981).

Determining observationally whether the mechanism we have described here operates in real systems will clearly be difficult. As the solar system example makes clear (Laskar & Robutel 1993), secular resonances in multiple-planet systems can be efficient sources of late-time obliquity (even Jupiter's low and presumably primordial obliquity is predicted to eventually increase; Saillenfest et al. 2020). Even when a system's current architecture is clearly nonresonant, excitation of obliquity could have occurred during a prior epoch of migration (Vokrouhlický & Nesvorný 2015). Young systems, observed during or shortly after the disk-embedded phase, provide the best opportunities. Our model suggests (and requires) that some circumplanetary disks around giant planets (Zhu 2015) ought to be found to be significantly tilted. In some cases the tilts may be large, approaching the critical misalignment for Kozai–Lidov oscillations (Martin et al. 2021). Planetary obliquities ought to be nonzero even for the youngest planets, and should be larger on average for planets at large orbital radii, which would have hosted more massive circumplanetary disks that are more efficient at communicating their misalignments to the planet. Lower-mass planets may also be favored, as their circumplanetary disks are expected to have a larger disk aspect ratio making them more unstable to tilting.

We thank Kaitlin Kratter and Zhaohuan Zhu for valuable discussions, and acknowledge support from NASA TCAN award 80NSSC19K0639.