The TESS–Keck Survey. IV. A Retrograde, Polar Orbit for the Ultra-low-density, Hot Super-Neptune WASP-107b

We used a Gaussian likelihood for the RV time series ( t , v _r ) given the model parameters Θ, and included a RV jitter term (σ_j ) to account for additional astrophysical or instrumental noise,

$$\begin{eqnarray}&&p({{\boldsymbol{v}}}_{r},{\boldsymbol{t}}| {\boldsymbol{\Theta }})=\displaystyle \prod _{i=1}^{N}\displaystyle \frac{1}{\sqrt{2\pi {\sigma }^{2}}}\exp \left[-\displaystyle \frac{{\left({v}_{r,i}-f({t}_{i},{\boldsymbol{\Theta }})\right)}^{2}}{2{\sigma }_{i}^{2}}\right],\end{eqnarray} \tag{ 1 }$$

where ${\sigma }_{i}^{2}={\sigma }_{\mathrm{RV},{\rm{i}}}^{2}+{\sigma }_{j}^{2}$. The model f(t_i , Θ) is given by

$$\begin{eqnarray}&&f({t}_{i},{\boldsymbol{\Theta }})=\mathrm{RM}({{\rm{t}}}_{i},{\boldsymbol{\theta }})+\gamma +\dot{\gamma }({{\rm{t}}}_{i}-{{\rm{t}}}_{0}),\end{eqnarray} \tag{ 2 }$$

where ${\boldsymbol{\Theta }}=({\boldsymbol{\theta }},\gamma ,\dot{\gamma })$ is the RM model parameters ( θ ) as well as an offset (γ) and slope ($\dot{\gamma }$) term which we added to approximate the reflex motion of the star and model any other systematic shift in RV throughout the transit (e.g., from noncrossed spots). The reference time t₀ is the time of the first observation (BJD).

RM(t_i , θ ) is the RM model described in Hirano et al. (2011). We assumed zero stellar differential rotation and adopted the transit parameters determined by Dai & Winn (2017), which came from a detailed analysis of K2 short-cadence photometry. We performed a simultaneous fit to the photometric and spectroscopic transit data using the same photometric data from K2 as in Dai & Winn (2017) to check for consistency. We obtained identical results for the transit parameters as they did, hence we opted to simply adopt their values, including their quadratic limb-darkening model. These transit parameters are all listed in Table 2. Our best-fit RV jitter is ${\sigma }_{j}={2.61}_{-0.51}^{+0.64}$ m s⁻¹, smaller than the jitter from the Keplerian fit to the full RV sample of ${3.9}_{-0.4}^{+0.5}$ m s⁻¹ (Piaulet et al. 2021). This is expected as the RM sequence covers a much shorter time baseline as compared to the full RV baseline, and as a result is only contaminated by short-term stellar noise sources such as granulation and convection.

Table 2. Adopted Parameters of the WASP-107 System

Parameter	Value	Unit	Source
Pb	5.7214742	days	1
tc	7584.329897 ± 0.000032	JD a	1
b	0.07 ± 0.07		1
iorb,b	${89.887}_{-0.097}^{+0.074}$	degrees	1
Rp /R⋆	0.14434 ± 0.00018		1
a/R⋆	18.164 ± 0.037		1
eb	0.06 ± 0.04		2
ωb	${40}_{-60}^{+40}$	degrees	2
Mb	30.5 ± 1.7	M⊕	2
Pc	${1088}_{-16}^{+15}$	days	2
ec	0.28 ± 0.07		2
ωc	$-{120}_{-20}^{+30}$	degrees	2
${M}_{c}\sin {i}_{\mathrm{orb},{\rm{c}}}$	0.36 ± 0.04	MJ	2
Teff	4245 ± 70	K	2
M*	${0.683}_{-0.016}^{+0.017}$	M⊙	2
R*	0.67 ± 0.02	R⊙	2
u1	0.6666 ± 0.0062		1
u2	0.0150 ± 0.0110		1

Note.

^a Days since JD 2,450,000. Sources: (1) Dai & Winn (2017); (2) Piaulet et al. (2020, submitted).

Download table as:  ASCIITypeset image

Table 3. WASP-107b Rossiter–McLaughlin Parameters

Parameter	MCMC CI	MAP value	Unit
Model Parameters
$\sqrt{v\sin {i}_{\star }}\cos \lambda$	$-{0.309}_{-0.154}^{+0.150}$	−0.30	a
$\sqrt{v\sin {i}_{\star }}\sin \lambda$	$-{0.126}_{-0.771}^{+0.808}$	−0.72	a
$\cos {i}_{s}$	$-{0.003}_{-0.681}^{+0.682}$	−0.56
γ	${0.80}_{-1.38}^{+1.36}$	0.97	m s−1
$\dot{\gamma }$	$-{20.83}_{-10.94}^{+11.05}$	−21.85	m s−2
σjit	${2.61}_{-0.51}^{+0.64}$	2.20	m s−1
$\mathrm{log}(| {v}_{{cb}}| )$	${0.89}_{-1.27}^{+1.18}$	2.17	a
Derived Parameters
∣λ∣	${118.1}_{-19.1}^{+37.8}$	112.63	degrees
$v\sin {i}_{\star }$	${0.45}_{-0.23}^{+0.72}$	0.61	km s−1
vcb	$-{7.74}_{-109.71}^{+7.33}$	−149.41	m s−1
i⋆	${28.17}_{-20.04}^{+40.38}$	7.06	degrees
∣ψ∣	${109.81}_{-13.64}^{+28.17}$	92.60	degrees

Note.

^a $v\sin {i}_{\star }$ is in km s⁻¹ and v_cb is in m s⁻¹.

Download table as:  ASCIITypeset image

The free parameters in the RM model are the sky-projected obliquity (λ), stellar inclination angle (i_⋆), and projected rotational velocity ($v\sin {i}_{\star }$). To the first order, the impact parameter b and sky-projected obliquity λ determine the shape of the RM signal, while $v\sin {i}_{\star }$ and R_p /R_⋆ set the amplitude. We adopted the parameterization ($\sqrt{v\sin {i}_{\star }}\cos \lambda$, $\sqrt{v\sin {i}_{\star }}\sin \lambda$) to improve the sampling efficiency and convergence of the Markov Chain Monte Carlo (MCMC). A higher order effect that becomes important when the RM amplitude is small is the convective blueshift, which we denote as v_cb (see Section 3.3 for more details). There are thus seven free parameters in our model: $\sqrt{v\sin {i}_{\star }}\cos \lambda$, $\sqrt{v\sin {i}_{\star }}\sin \lambda$, $\cos {i}_{\star }$, $\mathrm{log}(| {v}_{{cb}}| )$, γ, $\dot{\gamma }$, and σ_j . We placed a uniform hard-bounded prior on $v\sin {i}_{\star }\in [0,5]$ km s⁻¹ and on $\cos {i}_{\star }\in [0,1]$, and used a Jeffrey's prior for σ_j . All other parameters were assigned uniform priors.

The shape of the RM curve is also affected by processes on the surface of the star that broaden spectral lines, which affect the inferred RVs. In the Hirano et al. (2011) model, these processes are parameterized by γ_lw, the intrinsic line width, ζ, the line width due to macroturbulence, given by the Valenti & Fischer (2005) scaling relation

$$\begin{eqnarray}&&\zeta =\left(3.98+\displaystyle \frac{{T}_{\mathrm{eff}}-5770\,{\rm{K}}}{650\,{\rm{K}}}\right)\mathrm{km}\ {{\rm{s}}}^{-1},\end{eqnarray} \tag{ 3 }$$

and β, given by

$$\begin{eqnarray}&&\beta =\sqrt{\displaystyle \frac{2{k}_{B}{T}_{\mathrm{eff}}}{\mu }+{\xi }^{2}+{\beta }_{\mathrm{IP}}},\end{eqnarray} \tag{ 4 }$$

where ξ is the dispersion due to microturbulence and β_IP is the Gaussian dispersion due to the instrument profile, which we set to the HIRES line-spread function (LSF) (2.2 km s⁻¹). We tested having γ_lw, ξ, and ζ as free parameters in the model (with uniform priors) but only recovered the prior distributions for these parameters. Moreover we saw no change in the resulting posterior distribution for λ or $v\sin {i}_{\star }$. Because of this, we opted to instead adopt fixed nominal values of ξ = 0.7 km s⁻¹, γ_lw = 1 km s⁻¹, and ζ = 1.63 km s⁻¹ (from Equation (3) using T_eff from Table 2).

Convection in the stellar photosphere, caused by hotter bubbles of gas rising to the stellar surface and cooler gas sinking, results in a net blueshift across the stellar disk. This is because the rising (blueshifted) gas is hotter, and therefore brighter, than the cooler sinking (redshifted) gas. Since this net-blueshifted signal is directed at an angle normal to the stellar surface, the radial component seen by the observer is different in amplitude near the limb of the star compared to the center of the stellar disk, according to the stellar limb-darkening profile. Thus the magnitude of the convective blueshift blocked by the planet varies over the duration of the transit. The amplitude of this effect is ∼2 m s⁻¹, which is significant given the small amplitude of the RM signal we observe for WASP-107b (∼5.5 m s⁻¹).

For this reason we included the prescription of Shporer & Brown (2011) in the RM model, which is parameterized by the magnitude of the convective blueshift integrated over the stellar disk (v_cb ). This quantity is negative by convention. Since the possible value of v_cb could cover several orders of magnitude, we fit for $\mathrm{log}(| {v}_{{cb}}| )$ and set a uniform prior between −1 and 3. While we found that including v_cb has no effect on the recovered λ and $v\sin {i}_{\star }$ posteriors, we are able to rule out ∣v_cb ∣ > 450 m s⁻¹ at 99% confidence, and >250 m s⁻¹ at 95% confidence.

We first found the maximum a posteriori (MAP) solution by minimizing the negative log-posterior using Powells method (Powell 1964) as implemented in scipy.optimize.minimize (Virtanen et al. 2020). The MAP solution was then used to initialize an MCMC. We ran 8 parallel ensembles each consisting of 32 walkers for 10,000 steps using the python package emcee (Foreman-Mackey et al. 2013). We checked for convergence by requiring that both the Gelman–Rubin statistic (G–R; Gelman et al. 2003) was <1.001 across the ensembles (Ford 2006) and the autocorrelation time was <50 times the length of the chains (Foreman-Mackey et al. 2013).

The MAP values and central 68% confidence intervals (CI) computed from the MCMC chains are tabulated in Table 3, and the full posteriors for λ and $v\sin {i}_{\star }$ are shown in Figure 2. A prograde (∣λ∣ < 90°) orbit is ruled out at >99% confidence. An anti-aligned (135° < λ < 225°) orbit is allowed if $v\sin {i}_{\star }$ is small (0.26 ± 0.10 km s⁻¹), although a more polar aligned (but still retrograde) orbit with 90° < ∣λ∣ < 135° is more likely (if $v\sin {i}_{\star }\in [0.22,2.09]$ km s⁻¹, 90% CI). The true obliquity ψ will always be closer to a polar orientation than λ, since λ represents the minimum obliquity in the case where the star is viewed edge-on (i_⋆ = 90°). While an equatorial orbit that transits requires i_⋆ ∼ 90°, a polar orbit may be seen to transit for any stellar inclination.

To confirm that the signal we detected was not driven by correlated noise structures in the data, we performed a test using the cyclical residual permutation technique. We first calculated the residuals from the MAP fit to the original RV time series. We then shifted these residuals forward in time by one data point, wrapping at the boundaries, and added these new residuals back to the MAP model. This new "fake" data set was then fit again and the process was repeated N times where N = 22 is the number of data points in our RV time series. This technique preserves the red noise component, and permuting multiple times generates data sets that have the same temporal correlation but different realizations of the data. If we assume that the signal we detected is caused by a correlated noise structure, then we would expect to see the detected signal vanish or otherwise become significantly weaker across each permutation as that noise structure becomes asynchronous with the transit ephemeris. We found that the signal is robustly detected at all permutations, with and without including the convective blueshift (fixed to the original MAP value). The MAP estimate for λ tended to be closer to polar across the permutations compared as to the original fit, which is consistent with the posterior distribution estimated from the MCMC, but did not vary significantly. While this method is not appropriate for estimating parameter uncertainties (Cubillos et al. 2017), we conclude that our results are not qualitatively affected by correlated noise in our RV time series.

Spot-crossing events can also affect the RM curve since the planet would block a different amount of red/blueshifted light. Out of the nine transits observed by Dai & Winn (2017), a single spot-crossing event was seen in only three of the transits. Hence there is roughly a one in three chance that the transit we observed contained a spot-crossing event. As we did not obtain simultaneous high-cadence photometry, we do not know if or when such an event occurred. Judging from the durations (∼30 min) of the spot crossings observed by Dai & Winn (2017), this would only affect one or maybe two of our 15-minute exposures. While we dont see any significant outliers in our data set, these spots were only ∼10% changes on a ∼2% transit depth, amounting to an overall spot depth of ∼0.2%. Given our estimate of $v\sin {i}_{\star }\sim 0.5$ km s⁻¹, this suggests a spot-crossing event would produce a ∼1 m s⁻¹ RV anomaly, small compared to our measurement uncertainties (∼1.5 m s⁻¹) and the estimated stellar jitter (∼2.6 m s⁻¹). In other words, there is a roughly 33% chance that a spot-crossing event introduced an additional 0.5σ error on a single data point. If there were multiple spot-crossing events this anomaly would vary across the transit similar to other stellar-activity processes. In practice this introduces a correlated noise structure in the RV time series which our cyclical residual permutation test demonstrated is not significantly influencing our measurement of the obliquity or other model parameters. From this semi-analytic analysis we conclude that spot crossings are not a leading source of uncertainty in our model.

Given a constraint on $v\sin {i}_{\star }$ and v, we can constrain the stellar inclination i_⋆. Previous studies have found a range of estimates for the $v\sin {i}_{\star }$ of WASP-107. Anderson et al. (2017) found a value of 2.5 ± 0.8 km s⁻¹, whereas John Brewer (private communication) obtained a value of 1.5 ± 0.5 km s⁻¹ using the automated spectral synthesis modeling procedure described in Brewer et al. (2016). We note that the Specmatch-Emp (Yee et al. 2017) result for our HIRES spectrum only yields an upper bound for $v\sin {i}_{\star }$ of <2 km s⁻¹, as this technique is limited by the HIRES LSF. All three of these methods derive $v\sin {i}_{\star }$ by modeling the amount of line broadening present in the stellar spectrum, which in part comes from the stellar rotation. However these estimates may be biased from other sources of broadening which are not as well constrained in these models. Our RM analysis on the other hand incorporates a direct measurement of $v\sin {i}_{\star }$ by observing how much of the projected stellar rotational velocity is blocked by the transiting planet's shadow. Our RM analysis found $v\sin {i}_{\star }={0.45}_{-0.23}^{+0.72}$ km s⁻¹, lower than the spectroscopic estimates. We adopted this posterior for $v\sin {i}_{\star }$ to keep internal consistency.

The rotation period of WASP-107 has been estimated to be 17 ± 1 days from photometric modulations due to starspots rotating in and out of view (Anderson et al. 2017; Dai & Winn 2017; Močnik et al. 2017). We combined this rotation period with the stellar radius of 0.67 ± 0.02 R_⊙ inferred from the HIRES spectrum (Piaulet et al. 2021) using Specmatch-Emp (Yee et al. 2017) to constrain the tangential rotational velocity of v = 2π R_⋆/P_rot. We then used the statistically correct procedure described by Masuda & Winn (2020) and performed an MCMC sampling of v and $\cos {i}_{\star }$, using uniform priors for each, and using the posterior distribution for $v\sin {i}_{\star }$ obtained in the RM analysis as a constraint. Sampling both variables simultaneously correctly incorporates the nonindependence of v and $\cos {i}_{\star }$, since $v\leqslant v\sin {i}_{\star }$. We found that ${i}_{\star }={25.8}_{-15.4}^{+22.5}$ degrees (MAP value 7.1°), implying a viewing geometry of close to pole-on for the star. Thus any transiting configuration will necessarily imply a near-polar orbit, even for orbital solutions with λ near 180° (see Figure 3). It is worth mentioning that one of the three spot-crossing events observed by Dai & Winn (2017) occurred near the transit midpoint. This small stellar inclination implies that this spot must be at a relatively high latitude (90° − i_⋆) compared to that of our Sun, which has nearly all of its sunspots contained within ±30° latitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Knowledge of the stellar inclination i_⋆, the orbital inclination i_orb, and the sky-projected obliquity λ allows one to compute the true obliquity ψ, as these four angles are related by

$$\begin{eqnarray}&&\cos \psi =\cos {i}_{\mathrm{orb}}\cos {i}_{\star }+\sin {i}_{\mathrm{orb}}\sin {i}_{\star }\cos \lambda .\end{eqnarray} \tag{ 5 }$$

The resulting posterior distribution for the true obliquity ψ is shown in Figure 4. As expected, the true orbit is constrained to a more polar orientation than is implied by the wide posteriors on λ, due to the nearly pole-on viewing geometry of the star itself.

Zoom In Zoom Out Reset image size

Since the Hamiltonian ${ \mathcal H }$ does not depend on h_b , the quantity ${H}_{b}=\sqrt{1-{e}_{b}^{2}}\cos {i}_{b}$ is conserved. This leads to a periodic exchange of e_b and i_b , so long as the outer planet has an inclination greater than a critical value of ∼392 (Kozai 1962; Lidov 1962). These Kozai–Lidov cycles also require a slowly changing argument of perihelion, which may precess due to GR as is famously seen in the orbit of Mercury. This precession can suppress Kozai–Lidov cycles if fast enough, as is the case for HAT-P-11 and π Men (Xuan & Wyatt 2020; Yee et al. 2018). The precession rate from GR is given by

$$\begin{eqnarray}&&{\dot{\omega }}_{{GR}}=\displaystyle \frac{{{GM}}_{\star }}{{a}_{b}{c}^{2}}\displaystyle \frac{3{n}_{b}}{{G}_{b}^{2}},\end{eqnarray} \tag{ 8 }$$

which has an associated timescale of ${\tau }_{{GR}}=2\pi /\dot{\omega }\approx 42,500$ years for WASP-107b. The Kozai timescale (Kiseleva et al. 1998) is

$$\begin{eqnarray}&&{\tau }_{\mathrm{Kozai}}=\displaystyle \frac{2{P}_{c}^{2}}{3\pi {P}_{b}^{2}}\displaystyle \frac{{M}_{\star }}{{M}_{c}}{\left(1-{e}_{c}^{2}\right)}^{3/2}\approx 210,000\,\mathrm{yr},\end{eqnarray} \tag{ 9 }$$

five times longer. The condition for Kozai–Lidov cycles to be suppressed by relativistic precession is ${\tau }_{\mathrm{Kozai}}{\dot{\omega }}_{\mathrm{GR}}\gt 3$ (Fabrycky & Tremaine 2007), which the MAP minimum mass and orbital parameters WASP-107 c satisfy. This is nicely visualized in Figure 6 of Piaulet et al. (submitted), which shows the full posterior distributions of τ_Kozai and τ_GR. While the true mass of WASP-107 c is likely to be larger than the derived $M\sin {i}_{\mathrm{orb},{\rm{c}}}$, it would need to be ∼10 times larger for Kozai–Lidov oscillations to occur. This would imply a near face-on orbit of at most i_orb,c < 55. Such a face-on orbit is unlikely but is still plausible if it is aligned with the rotation axis of the star, given our constraints on the stellar inclination angle in Section 3.5.

An alternative explanation for the high obliquity of WASP-107b is nodal precession, as was proposed for HAT-P-11b (Yee et al. 2018) and for π Men c (Xuan & Wyatt 2020). In this scenario the outer planet must have an obliquity greater than half that of the inner planet, which in this case would require ψ_c ∼ 55°. Then the longitude of ascending node Ω_b evolves in a secular manner according to Yee et al. (2018),

$$\begin{eqnarray}&&\displaystyle \frac{d{{\rm{\Omega }}}_{b}}{{dt}}=\displaystyle \frac{\partial { \mathcal H }}{\partial {H}_{b}}=\displaystyle \frac{{n}_{b}}{8}\displaystyle \frac{{M}_{c}}{{M}_{\star }}{\left(\displaystyle \frac{{a}_{b}}{{a}_{c}\sqrt{1-{e}_{c}^{2}}}\right)}^{3}\left(\displaystyle \frac{15-9{G}_{b}^{2}}{{G}_{b}^{2}}\right){H}_{b}.\end{eqnarray} \tag{ 10 }$$

The associated timescale ${\tau }_{{{\rm{\Omega }}}_{b}}=2\pi /\dot{{{\rm{\Omega }}}_{b}}$ is only about 2 Myr, much shorter than the age of the system. Yee et al. (2018) pointed out that such a precession will cause the relative inclination of the two planets to oscillate between ≈ ψ_c ± ψ_c . Thus at certain times the observer may see a highly misaligned orbit (ψ_b ∼ 2ψ_c ) for the inner planet, while at other times the observer may see an aligned orbit (ψ_b = 0).

We examined this effect by running a 3D N-body simulation in REBOUND (Rein & Liu 2012). We initialized planet c with an obliquity of 60° (which sets the maximum obliquity planet b can obtain, ∼2ψ_c = 120°) and planet b with an obliquity of 0° (aligned, prograde orbit). We included the effects of GR and tides using the gr and modify_orbits_forces features of REBOUNDx (Kostov et al. 2016; Tamayo et al. 2019) and used the the WHFast integrator (Rein & Tamayo 2015) to evolve the system forward in time for 10 Myr.

Figure 5 shows that over these 10 Myr ψ_b oscillates in the range of 0°–120° due to the precession of Ω_b . Thus nodal precession can easily produce high relative inclinations, despite Kozai–Lidov oscillations being suppressed by GR. A configuration like what is observed today in which the inner planet is misaligned on a polar, yet slightly retrograde orbit is attainable at times during this cycle where the mutual inclination is at or near its maximum. The obliquity is ≳80% the amplitude from nodal precession (∼2ψ_c ) approximately one-third of the time (bottom panel in Figure 6). Therefore, even though the observed obliquity depends on when during the nodal precession cycle the system is observed, there is a decent chance of observing ψ_b near its maximum.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the simulation we ran, WASP-107b is only seen by an observer to be in a transiting geometry about 2.8% of the time. Xuan & Wyatt (2020) did a more detailed calculating accounting for the measured mutual inclination and found that the dynamical transit probability for π Men c and HAT-P-11b is of order 10%–20%. However, as Xuan & Wyatt (2020) point out, this does not affect the population-level transit likelihood since the overall orientations of extrasolar systems can still be treated as isotropic. It merely suggests that a system with a transiting distant giant planet may be harboring a nodally precessing inner planet that just currently happens to be nontransiting.

Both Kozai–Lidov and nodal precession require a large mutual inclination in order for the inner planet to reach polar orientations. The origin of this large mutual inclination may be hidden in the planet's formation history, or perhaps was caused by a planet–planet scattering event with an additional companion that was ejected from the system. This could also explain the moderately eccentric orbit of WASP-107 c (Piaulet et al. 2021). Indeed a significant mutual inclination is observed for the inner and outer planets of the HAT-P-11 and π Men systems (Xuan & Wyatt 2020), although the inner planet in π Men is only slightly misaligned with λ = 24 ± 4.1 degrees (Kunovac Hodžić et al. 2021), while HAT-P-11b has $\lambda ={103}_{-10}^{+26}$ degrees (Winn et al. 2010b).

As more close-in Neptunes with distant giant companions are discovered, the distribution of observed obliquities for the inner planet will help determine if we are indeed simply seeing many systems undergoing nodal precession but at different times during the precession cycle. If so, we might observe a sky-projected obliquity distribution that resembles the bottom panel of Figure 6. However, we may instead be observing two classes of close-in Neptunes: ones aligned with their host stars and ones in polar or near-polar orbits (see the top panel of Figure 6). This suggests an alternative mechanism that favors either polar orbits or aligned orbits depending on the system architecture.

Recently, Petrovich et al. (2020) showed that, even for ψ_c ∼ 0°, a resonance encountered as the young protoplanetary disk dissipates can excite an inner planet to high obliquities, even favoring a polar orbit given appropriate initial conditions. To summarize the model, consider a system with a close-in planet and a distant (few astronomical units) giant planet, like WASP-107, after the disk interior to the outer planet has been cleared but the disk exterior remains. The external gaseous disk induces a nodal precession of the outer planet at a rate proportional to the disk mass (Equation (10) with b ↦ c and c ↦ disk). The outer planet still induces a nodal precession on the inner planet according to Equation (10). If at first the rate dΩ_c /dt > dΩ_b /dt, then as the disk dissipates (and M_disk decreases) the precession rate for planet c will decrease until it matches the precession rate of the inner planet. At this point the system will pass through a secular resonance, driving an instability which tilts the inner planet to a high obliquity; a small initial obliquity of a few degrees can quickly reach 90°. Additionally, depending on the relative strength of the stellar quadrupole moment and GR effects, the inner planet may obtain a high eccentricity (if GR is unimportant), a modest eccentricity (if GR is important), or a circular orbit (if GR dominates). Tidal forces can circularize the orbit, although the planet may retain a detectable eccentricity even after several gigayears. This process well explains the polar, close-in, and eccentric orbits of small planets like HAT-P-11b. Nodal precession alone is unable to explain the eccentricity of such planets.

Given the planet and stellar properties of the WASP-107 system, we calculated the instability criteria developed in Petrovich et al. (2020). The steady-state evolution of the system can be inferred by comparing the relative strength of GR (η_GR) with the stellar quadrupole moment (η_⋆). We found that η_GR > η_⋆ + 6 at 99.76% confidence, η_⋆ + 6 > η_GR > 4 at 0.155% confidence, and η_GR < 4 at 0.084% confidence (i.e., η_GR ∼ 30–80 and η_⋆ ∼ 1). Thus WASP-107b is stable against eccentricity instabilities and lives in the polar, circular region of parameter space in Figure 4 of Petrovich et al. (2020).

We calculated the final obliquity of WASP-107b using the procedure outlined in Petrovich et al. (2020), incorporating the uncertainties in $M\sin {i}_{\mathrm{orb},{\rm{c}}}$ and P_c and integrating over all possible initial obliquities for the outer planet. Evaluating their Equation (3), we found that the resonance that drives the inner planet to high obliquities is always crossed. We calculated the adiabatic parameter x_ad ≡ τ_disk/τ_adia from the disk dispersal timescale and the adiabatic time (their Equation (7)), taking τ_disk to be 1 Myr. In the orbital configurations where x_ad > 1 (adiabatic crossing) we computed the final obliquity from their Equation (12) (I_crit). Otherwise, the final obliquity was set to I_non−ad from their Equation (15).

The resulting probability of the final obliquity of WASP-107b is 7.6% for a nonpolar (but oblique) orbit and 92.4% for a polar orbit. A polar orbit is likely if the outer planet's orbit is inclined at least ∼8°, and is guaranteed for ψ_init,c ≳ 25°. In an equivalent parameterization, Petrovich et al. (2020) explicitly predict a polar orbit for WASP-107b if the mass and semiminor axis of WASP-107 c satisfy ${\left({b}_{c}/2\mathrm{au}\right)}^{3}\gt ({M}_{c}/0.5\,{M}_{J})$. Since we only have a constraint on $M\sin {i}_{\mathrm{orb},{\rm{c}}}$, this condition is satisfied if i_orb,c ∈ [60°–90°]. Such a viewing geometry, in conjunction with an obliquity of ψ_c > 25°, is plausible given the likely stellar orientation (Section 3.5).

A key deviation from this model is that while the orbit of WASP-107b is indeed close to polar, it is quite definitively retrograde. In the disk dispersal-driven tilting scenario, the inner planet approaches a ψ = 90° polar orbit from below and stops at ψ_b = 90°. In order to reach a super-polar/retrograde orbit, WASP-107 c must have a significant obliquity, either primordial from formation or through a scattering event (Petrovich et al. 2020). As we alluded to in Section 4.2, a scattering event could also explain the moderate eccentricities of the outer giants WASP-107 c and HAT-P-11 c, and could easily give WASP-107 c a high enough obliquity to guarantee a polar/super-polar configuration for WASP-107b (Huang et al. 2017). In fact a scattering event is more likley to produce the modest obliquity for planet c needed to produce a super-polar orbit under the disk dispersal framework than it is to produce the large (ψ_c ≳ 40–50°) obliquity needed to excite either Kozai–Lidov or nodal precession cycles.