The TESS-Keck Survey IV: A Retrograde, Polar Orbit for the Ultra-Low-Density, Hot Super-Neptune WASP-107b

We measured the Rossiter-McLaughlin effect of WASP-107b during a single transit with Keck/HIRES. We found the sky-projected inclination of WASP-107b's orbit, relative to its host star's rotation axis, to be $|\lambda| = {118}^{+38}_{-19}$ degrees. This confirms the misaligned/polar orbit that was previously suggested from spot-crossing events and adds WASP-107b to the growing population of hot Neptunes in polar orbits around cool stars. WASP-107b is also the fourth such planet to have a known distant planetary companion. We examined several dynamical pathways by which this companion could have induced such an obliquity in WASP-107b. We find that nodal precession and disk dispersal-driven tilting can both explain the current orbital geometry while Kozai-Lidov cycles are suppressed by general relativity. While each hypothesis requires a mutual inclination between the two planets, nodal precession requires a much larger angle which for WASP-107 is on the threshold of detectability with future Gaia astrometric data. As nodal precession has no stellar type dependence, but disk dispersal-driven tilting does, distinguishing between these two models is best done on the population level. Finding and characterizing more extrasolar systems like WASP-107 will additionally help distinguish whether the distribution of hot-Neptune obliquities is a dichotomy of aligned and polar orbits or if we are uniformly sampling obliquities during nodal precession cycles.


INTRODUCTION
WASP-107b is a close-in (P = 5.72 days) super-Neptune orbiting the cool K-dwarf WASP-107. Originally discovered via the transit method by WASP-South, WASP-107b was later observed by K2 in Campaign 10 (Howell et al. 2014). These transits revealed a radius close to that of Jupiter, R b = 10.8 ± 0.34 R ⊕ = * NSF Graduate Research Fellow 0.96 ± 0.03 R J (Dai & Winn 2017;Močnik et al. 2017;Piaulet et al. 2021). However, follow-up radial velocity (RV) measurements with the CORALIE spectrograph demonstrated a mass of just 38 ± 3 M ⊕ , meaning this Jupiter-sized planet has just onetenth its density. Higher-precision RVs from Keck/High Resolution Echelle Spectrometer (HIRES) suggested an even lower mass of 30.5 ± 1.7 M ⊕ (Piaulet et al. 2021). This low density challenges the standard core-accretion model of planet formation. If runaway accretion brought WASP-107b to a gas-to-core mass ratio of ∼ 3 but was stopped prematurely before growing to gas giant size, orbital dynamics and/or migration may have played a significant role in this system (Piaulet et al. 2021). Alternatively WASP-107b's radius may be inflated from tidal heating, which would allow a lower gas-to-core ratio consistent with core accretion (Millholland et al. 2020).
With a low density, large radius, and hot equilibrium temperature, WASP-107b's large atmospheric scale height makes it a prime target for atmospheric studies. Indeed analyses of transmission spectra obtained with the Hubble Space Telescope (HST)/WFC3 have detected water amongst a methane-depleted atmosphere (Kreidberg et al. 2018). WASP-107b was the first exoplanet to be observed transiting with excess absorption at 10830 Å, an absorption line of a metastable state of neutral helium indicative of an escaping atmosphere (Oklopčić & Hirata 2018). These observations suggest that WASP-107b's atmosphere is photoevaporating at a rate of a few percent in mass per billion years (Spake et al. 2018;Allart et al. 2019;Kirk et al. 2020).
The orbit of WASP-107b is suspected to be misaligned with the rotation axis of its host star. The angle between the star's rotation axis and the normal to the planet's orbital plane, called the stellar obliquity ψ (or just obliquity), was previously constrained by observations of WASP-107b passing over starspots as it transited (Dai & Winn 2017). As starspots are regions of reduced intensity on the stellar photosphere that rotate with the star, this is seen as a bump of increased brightness in the transit light curve. By measuring the time between spot-crossing events across successive transits, combined with the absence of repeated spot crossings, Dai & Winn (2017) were able to constrain the sky-projected obliquity, λ, of WASP-107b to λ ∈[40-140] deg. Intriguingly, long-baseline RV monitoring of the system with Keck/HIRES has revealed a distant (P c ∼ 1100 days) massive (M sin i orb,c = 115 ± 13 M ⊕ ) planetary companion, which may be responsible for this present day misaligned orbit through its gravitational influence on WASP-107b (Piaulet et al. 2021).
The sky-projected obliquity can also be measured spectroscopically. The Rossiter-McLaughlin (RM) effect refers to the anomalous Doppler-shift caused by a transiting planet blocking the projected rotational velocities across the stellar disk (McLaughlin 1924;Rossiter 1924). If the planet's orbit is aligned with the rotation of the star (prograde), its transit will cause an anomalous redshift followed by an anomalous blueshift. A antialigned (retrograde) orbit will cause the opposite to occur.
Following the first obliquity measurement by Queloz et al. (2000), the field saw measurements of 10 exoplanet obliquities over the next 8 years that were all consistent with aligned, prograde orbits. After a few misaligned systems had been discovered (e.g., Hébrard et al. 2008), a pattern emerged with hot Jupiters on highly misaligned orbits around stars hotter than about 6250 K (Winn et al. 2010a). This pattern elicited several hypotheses such as damping of inclination by the convective envelope of cooler stars (Winn et al. 2010a) or magnetic realignment of orbits during the T Tauri phase (Spalding & Batygin 2015).
In this paper we present a determination of the obliquity of WASP-107b from observations of the RM effect (Section 2). These observations were acquired under the TESS-Keck Survey (TKS), a collaboration between scientists at the University of California, the California Institute of Technology, the University of Hawai'i, and NASA. TKS is organized through the California Planet Search with the goal of acquiring substantial RV follow-up observations of planetary systems discovered by TESS (Dalba et al. 2020). TESS observed four transits of WASP-107b (TOI 1905) in Sector 10. An additional science goal of TKS is to measure the obliquities of interesting TESS systems. WASP-107b, which is already expected to have a significant obliquity (Dai & Winn 2017), is an excellent target for an RM measurement with HIRES.
In Section 3 we confirm a misaligned orientation; in fact, we found a polar/retrograde orbit. This adds WASP-107b to the growing population of hot Neptunes in polar orbits around cool stars. We explored possible mechanisms that could be responsible for this misalignment in Section 4. Lastly in Section 5 we summarized our findings and discussed the future work needed to better understand the obliquity distribution for small planets around cool stars. We observed the RM effect for WASP-107b during a transit on 2020 February 26 (UTC) with HIRES (Vogt et al. 1994) on the Keck I Telescope on Maunakea. Our HIRES observations covered the full transit duration (∼ 2.7 hr) with a ∼ 1 hour baseline on either side. We used the "C2" decker (14 × 0. 861, R = 45, 000) and integrated until the exposure meter reached 60,000 counts (signal-to-noise ratio (S/N) ∼ 100 per reduced pixel, 15 minutes) or readout after 15 minutes. The spectra were reduced using the standard procedures of the California Planet Search (Howard et al. 2010), with the iodine cell serving as the wavelength reference (Butler et al. 1996). In total we obtained 22 RVs, 12 of which were in transit (Table 1).
Visually inspecting the observations (Fig. 1) shows an anomalous blueshift following the transit ingress, followed by an anomalous redshift after the transit mid-  point, 1 , indicating a retrograde orbit. The asymmetry and low-amplitude of the signal constrain the orientation to a near-polar alignment, but whether the orbit is polar or anti-aligned is somewhat degenerate with the value of v sin i . The expected RM amplitude is v sin i (R p /R ) 2 ∼ 40 m s −1 , using previous estimates of R p /R = 0.144 (Dai & Winn 2017) and v sin i ∼ 2 km s −1 (e.g., Anderson et al. 2017). The signal we detected with HIRES is only ∼ 5.5 m s −1 in amplitude. Dai & Winn (2017) found the transit impact parameter to be nearly zero, therefore the small RM amplitude suggests either a much lower v sin i than was spectroscopically inferred (see Section 3.5), a near-polar orbit, or both.

Rossiter-McLaughlin Model
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 astrophys-ical or instrumental noise, is given by where Θ = (θ, γ,γ) is the RM model parameters (θ) as well as an offset (γ) and slope (γ) 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 0 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 σ j = 2.61 +0.64 −0.51 m s −1 , smaller than the jitter from the Keplerian fit to the full RV sample of 3.9 +0.5 −0.4 m s −1 (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.
The free parameters in the RM model are the skyprojected obliquity (λ), stellar inclination angle (i ), and projected rotational velocity (v sin i ). To first order, the impact parameter b and sky-projected obliquity λ determine the shape of the RM signal, while v sin i and R p /R set the amplitude. We adopted the param- 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 v cb (see Section 3.3 for more details). There are thus seven free parameters in our model: √ v sin i cos λ, √ v sin i sin λ, cos i , log(|v cb |), γ,γ, and σ j . We placed a uniform hard-bounded prior on v sin i ∈ [0, 5] km s −1 and on cos i ∈ [0, 1], and used a Jeffrey's prior for σ j . All other parameters were assigned uniform priors.

Micro/Macroturbulence Parameters
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 and β, given by 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 −1 ). 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 . Because of this, we opted to instead adopt fixed nominal values of ξ = 0.7 km s −1 , γ lw = 1 km s −1 , and ζ = 1.63 km s −1 (from Eq. 3 using T eff from Table 2).

Convective blueshift
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 limbdarkening 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 −1 , which is significant given the small amplitude of the RM signal we observe for WASP-107b (∼ 5.5 m s −1 ).
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 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 posteriors, we are able to rule out |v cb | > 450 m s −1 at 99% confidence, and > 250 m s −1 at 95% confidence.

degrees
a v sin i is in km s −1 and v cb is in m s −1 .

Evidence for a Retrograde/Polar Orbit
We first found the maximum a posteriori (MAP) solution by minimizing the negative log-posterior using Powell's 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 are shown in Fig. 2 .09] km s −1 , 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" dataset 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 datasets 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 don't see any significant outliers in our dataset, 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 ∼ 0.5 km s −1 this suggests a spot-crossing event would produce a ∼ 1 m s −1 RV anomaly, small compared to our measurement uncertainties (∼ 1.5 m s −1 ) and the estimated stellar jitter (∼ 2.6 m s −1 ). 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.

Constraints on the Stellar Inclination
Given a constraint on v sin i and v, we can constrain the stellar inclination i . Previous studies have found a range of estimates for the v sin i of WASP-107. Anderson et al. (2017) found a value of 2.5 ± 0.8 km s −1 , whereas John Brewer (private communication) obtained a value of 1.5 ± 0.5 km s −1 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 of < 2 km s −1 , as this technique is limited by the HIRES PSF. All three of these methods derive v sin i 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 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 = 0.45 +0.72 −0.23 km s −1 , lower than the spectroscopic estimates. We adopted this posterior for v sin i 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 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 v = 2πR /P rot . We then used the statistically correct procedure described by  and performed an MCMC sampling of v and cos i , using uniform priors for each, and using the posterior distribution for v sin i obtained in the RM analysis as a constraint. Sampling both variables simultaneously correctly incorporates the nonindependence of v and cos i , since v ≤ v sin i . We found that i = 25.8 +22.5 −15.4 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 Fig. 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.
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 cos ψ = cos i orb cos i + sin i orb sin i cos λ.
The resulting posterior distribution for the true obliquity ψ is shown in Fig. 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.

DYNAMICAL HISTORY
How did WASP-107b end up in a slightly retrograde, nearly polar orbit? To explore this question, we examined the orbital dynamics of the WASP-107 system considering the new discovery of a distant, giant companion WASP-107c (Piaulet et al. 2021). As in Mardling (2010), Yee et al. (2018), and Xuan & Wyatt (2020), we can understand the evolution of the WASP-107 system by examining the secular three-body Hamiltonian. Assuming the inner planet is a test particle (i.e., M b √ a b M c √ a c ), and since a b /a c 1, we can approximate the Hamiltonian by expanding to quadrupole order in semimajor axis ratio where the last term is the addition from general relativity (GR) and n b = 2π/P b . The quantities G and H are the canonical Delaunay variables where the double-arrow (↔) symbolizes conjugate variables, ω b is the argument of perihelion of the inner planet, Ω b is the longitude of ascending node of the inner planet, and i b is the inclination of the inner planet with respect to the invariant plane. The invariant plane is the plane normal to the total angular momentum bmtor, which to good approximation is simply the orbital plane of the outer planet (since angular momentum is ∝ M a 1/2 ). With this approximation, i b is the relative inclination between the two planets.

Kozai-Lidov oscillations
Since the Hamiltonian H does not depend on h b , the quantity H b = 1 − e 2 b 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 ∼ 39. • 2 (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 bẏ which has an associated timescale of τ GR = 2π/ω ≈ 42, 500 years for WASP-107b. The Kozai timescale (Kiseleva et al. 1998) is five times longer. The condition for Kozai-Lidov cycles to be suppressed by relativistic precession is τ KozaiωGR > 3 (Fabrycky & Tremaine 2007), which the MAP minimum mass and orbital parameters WASP-107c 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-107c is likely to be larger than the derived M sin i orb,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 < 5. • 5. Such a faceon 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.

Nodal precession
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), (10) The associated timescale τ Ω b = 2π/Ω 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. Fig. 5 shows that over these 10 Myr ψ b oscillates in the range 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 Fig. 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.
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-107c (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 (Hodžić et al. submitted), while HAT-P-11b has λ = 103 +26 −10 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 Fig. 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 Fig. 6). This suggests an alternative mechanism that favors either polar orbits or aligned orbits depending on the system architecture.

Disk dispersal-driven tilting
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 (Eq. 10 with b → c and c → disk). The outer planet still induces a nodal precession on the inner planet according to Eq. 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  sky-projected obliquity as the azimuthal coordinate and normalized orbital distance as the radial coordinate, for <100 M⊕ planets around stars with T eff < 6250 K (similar mass planets around hotter stars are shown as faded gray points). The red point is WASP-107b. Other noteworthy systems are shown with various colors and markers (see Section 1 for references). Data compiled from TEPCat as of 2020 October (Southworth 2011). Only WASP-107, HAT-P-11, and π Men have distant giant companions detected. Kepler-56 (Huber et al. 2013) is another similar system but is not included in this plot as it is an evolved massive star. Bottom: the fraction of a nodal precession cycle spent in a given obliquity bin (left). The true obliquity ψ is assumed to vary as cos[(π/2)ψ(t)/ψmax] = sin 2 (πt/τ ), where t ∈ [0, τ = 1]. This recreates the shape of the oscillating inclination in Fig. 5. The amplitude ψmax is twice the outer planet's inclination which is plotted for three different distributions (shown on the right): uniform between [0 • , 90 • ] (gray), uniform between [40 • , 60 • ] (red), and using the von-Mises Fisher distribution from  calculated in a hierarchical manner incorporating their posterior distribution for the shape parameter σ for all. In all three cases the true obliquity is shown as a dashed histogram. The sky-projected obliquity is computed given a transiting geometry (i orb,b = 90 • ) and is marginalized over stellar inclination angle (solid histogram). Mp < 100 M⊕ planets with observed sky-projected obliquities are shown as a filled histogram for comparison. Note that while the gray and black predictions are relatively similar, an excess of polar orbits can be observed if the mutual inclination distribution is clustered around ∼ 40-60 • . 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 Fig. 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 orb,c and P c and integrating over all possible initial obliquities for the outer planet. Evaluating their Eq. (3), we found that the res-onance 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 Eq. 7), taking τ disk to be 1 Myr. In the orbital configurations where x ad > 1 (adiabatic crossing) we computed the final obliquity from their Eq. (12) (I crit ). Otherwise, the final obliquity was set to I non-ad from their Eq. (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-107c satisfy (b c /2 AU) 3 > (M c /0.5 M J ). Since we only have a constraint on M sin i orb,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-107c 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-107c and HAT-P-11c, and could easily give WASP-107c 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.

DISCUSSION AND CONCLUSION
We observed the RM effect during a transit of WASP-107b on 2020 February 26, from which we derived a near-polar and retrograde orbit as well as a low stellar v sin i . This low v sin i implies that we are viewing the star close to one of its poles, reinforcing the nearpolar orbital configuration of WASP-107b. However, we are unable to conclusively say how WASP-107b acquired such an orbit. Nodal precession or disk dispersal-driven tilting are both plausible mechanisms for producing a polar orbit, while Kozai-Lidov oscillations may be possible but only for a very narrow range of face-on or-bital geometries for WASP-107c. RV observations (Piaulet et al. 2021) as well as constraints on the velocity of the escaping atmosphere of WASP-107b (e.g., Allart et al. 2019, Kirk et al. 2020, Spake, J. J. et al. 2020 are consistent with a circular orbit. The eccentricity damping timescale due to tidal forces is only ∼ 60 Myr (Piaulet et al. 2021), so this is not unexpected. While a circular orbit does not rule out any of these pathways, only disk dispersal-driven tilting can explain both the eccentric and polar orbit of WASP-107b's doppelganger HAT-P-11 b.
Since all three scenarios depend on the obliquity of the outer giant planet, measuring the mutual inclination of planet b and c is essential to understand the dynamics of this system. This has been done for similar system architectures such as HAT-P-11 (Xuan & Wyatt 2020) and π Men (Xuan & Wyatt 2020;De Rosa et al. 2020) by observing perturbations in the astrometric motion of the star due to the gravitational tugging of the distant giant planet, using data from Hipparcos and Gaia. Unfortunately WASP-107 is significantly fainter (V = 11.5; Anderson et al. 2017) and barely made the cutoff in the Tycho-2 catalog of Hipparcos (90% complete at V=11.5; Høg et al. 2000). The poor Hipparcos astrometric precision, combined with the small angular scale of the orbit of WASP-107 on the sky (10 -30 µas), prevents a detection of the outer planet using astrometry. Assuming future Gaia data releases have the same astrometric precision as in DR2 (44 µas for WASP-107), WASP-107c will be at the threshold of detectability using the full five-year astrometric time series.
On the population level, the disk dispersal-driven model favors low-mass and slowly rotating stars due to its dependence on the stellar quadrupole moment, and also can explain eccentric polar orbits. Since nodal precession has no stellar type preference nor a means of exciting eccentric orbits, measuring the obliquities and eccentricities for a population of close-in Neptunes will be essential for distinguishing which process is the dominant pathway to polar orbits. Additionally a large population is needed to determine if the overall distribution of planet obliquities is consistent with catching systems at different stages of nodal precession, or if there are indeed two distinct populations of aligned or polar close-in Neptunes. As these models all depend on the presence of an outer giant planet, long-baseline RV surveys will be instrumental for discovering the nature of any perturbing companions (e.g. Rosenthal et al. submitted). Moreover RV monitoring of systems with small planets that already have measured obliquities, but do not have mass constraints or detected outer companions, will further expand this population. Recent examples of such systems include Kepler-408b (Kamiaka et al. 2019), AU Mic b (Palle et al. 2020), HD 63433 (b, Mann et al. 2020and c, Dai et al. 2020), K2-25b (Stefánsson et al. 2020), and DS Tuc b (Montet et al. 2020;Zhou et al. 2020). Comparing the proportions of systems with and without companions which have inner aligned or misaligned planets will further illuminate the likelihood of these different dynamical scenarios.