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

Ryan A. Rubenzahl NSF Graduate Research Fellow  Department of Astronomy,
California Institute of Technology, Pasadena, CA 91125, USA  Fei Dai Division
of Geological and Planetary Sciences, California Institute of Technology,
Pasadena, CA 91125, USA  Andrew W. Howard Department of Astronomy, California
Institute of Technology, Pasadena, CA 91125, USA  Ashley Chontos NSF Graduate
Research Fellow  Institute for Astronomy, University of Hawai‘i, 2680 Woodlawn
Drive, Honolulu, HI 96822, USA  Steven Giacalone Department of Astronomy,
University of California Berkeley, Berkeley, CA 94720, USA  Jack Lubin
Department of Physics & Astronomy, University of California Irvine, Irvine, CA
92697, USA  Lee J. Rosenthal Department of Astronomy, California Institute of
Technology, Pasadena, CA 91125, USA  Howard Isaacson Department of Astronomy,
University of California Berkeley, Berkeley CA 94720, USA  Centre for
Astrophysics, University of Southern Queensland, Toowoomba, QLD, Australia
Natalie M. Batalha Department of Astronomy and Astrophysics, University of
California, Santa Cruz, CA 95060, USA  Ian J. M. Crossfield  Department of
Physics & Astronomy, University of Kansas, 1082 Malott, 1251 Wescoe Hall Dr.,
Lawrence, KS 66045, USA  Courtney Dressing Department of Astronomy, University
of California Berkeley, Berkeley CA 94720, USA  Benjamin Fulton NASA Exoplanet
Science Institute/Caltech-IPAC, MC 314-6, 1200 E. California Blvd., Pasadena,
CA 91125, USA  Daniel Huber Institute for Astronomy, University of Hawai‘i,
2680 Woodlawn Drive, Honolulu, HI 96822, USA  Stephen R. Kane Department of
Earth and Planetary Sciences, University of California, Riverside, CA 92521,
USA  Erik A Petigura Department of Physics & Astronomy, University of
California Los Angeles, Los Angeles, CA 90095, USA  Paul Robertson Department
of Physics & Astronomy, University of California Irvine, Irvine, CA 92697, USA
Arpita Roy Space Telescope Science Institute, 3700 San Martin Dr., Baltimore,
MD 21218, USA  Department of Physics and Astronomy, Johns Hopkins University,
3400 N Charles St, Baltimore, MD 21218, USA  Lauren M. Weiss Institute for
Astronomy, University of Hawai‘i, 2680 Woodlawn Drive, Honolulu, HI 96822, USA
Corey Beard Department of Physics & Astronomy, University of California
Irvine, Irvine, CA 92697, USA  Michelle L. Hill  Department of Earth and
Planetary Sciences, University of California, Riverside, CA 92521, USA  Andrew
Mayo Department of Astronomy, University of California Berkeley, Berkeley CA
94720, USA  Teo Močnik Gemini Observatory/NSF’s NOIRLab, 670 N. A’ohoku Place,
Hilo, HI 96720, USA  Joseph M. Akana Murphy NSF Graduate Research Fellow
Department of Astronomy and Astrophysics, University of California, Santa
Cruz, CA 95064, USA  Nicholas Scarsdale Department of Astronomy and
Astrophysics, University of California, Santa Cruz, CA 95060, USA

###### Abstract

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.

††software:  emcee (Foreman-Mackey et al., 2013), corner.py (Foreman-Mackey,
2016), REBOUND/REBOUND/x (Rein & Liu, 2012; Tamayo et al., 2019), SciPy
(Virtanen et al., 2020), NumPy (Harris et al., 2020), Matplotlib (Hunter,
2007)

##  1 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\pm
0.34~{}R_{\oplus}=0.96\pm 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\pm
3~{}M_{\oplus}$ (Anderson et al., 2017), meaning this Jupiter-sized planet has
just one-tenth its density. Higher-precision RVs from Keck/High Resolution
Echelle Spectrometer (HIRES) suggested an even lower mass of $30.5\pm
1.7~{}M_{\oplus}$ (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 $\sim 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 $\psi$ (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, $\lambda$, of WASP-107b to
$\lambda\in$[40–140] deg. Intriguingly, long-baseline RV monitoring of the
system with Keck/HIRES has revealed a distant ($P_{c}\sim 1100$ days) massive
($M\sin i_{\text{orb},c}=115\pm 13~{}M_{\oplus}$) 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 anti-aligned
(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).

More recently a number of exoplanets have been found on misaligned orbits
around cooler stars, such as the hot Jupiter WASP-8b (Queloz et al., 2010;
Bourrier et al., 2017), as well as lower-mass hot Neptunes like HAT-P-11b
(Winn et al., 2010b), Kepler-63b (Sanchis-Ojeda et al., 2013), HAT-P-18b
(Esposito, M. et al., 2014), GJ 436b (Bourrier et al., 2018), and HD 3167 c
(Dalal et al., 2019). Strikingly, all of these exoplanets are on or near polar
orbits. Some of these systems have recently had distant, giant companions
detected (e.g. HAT-P-11c; Yee et al., 2018), hinting that these obliquities
arise from multibody planet-planet dynamics.

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.

Table 1: Radial Velocities of WASP-107

Time | RV | $\sigma_{\rm{RV}}$ | Exposure time
---|---|---|---
(BJD${}_{\text{TDB}}$) | (m s-1) | (m s-1) | (s)
2458905.90111 | 5.05 | 1.50 | 900
2458905.91189 | 6.43 | 1.42 | 883
2458905.92247 | 0.14 | 1.49 | 862
2458905.93288 | -1.35 | 1.65 | 844
2458905.94266 | -0.25 | 1.45 | 783
2458905.95204 | -5.28 | 1.44 | 745
2458905.96141 | -2.40 | 1.37 | 797
2458905.97098 | -3.40 | 1.46 | 754
2458905.98004 | 2.45 | 1.37 | 727
2458905.98927 | -5.52 | 1.45 | 780
2458905.99888 | 2.07 | 1.48 | 792
2458906.00848 | 4.21 | 1.37 | 776
2458906.01796 | -0.58 | 1.38 | 775
2458906.02768 | 0.83 | 1.47 | 817
2458906.03780 | 3.07 | 1.49 | 836
2458906.04780 | -3.01 | 1.26 | 818
2458906.05771 | 0.02 | 1.45 | 796
2458906.06752 | -3.72 | 1.49 | 795
2458906.07703 | 3.61 | 1.33 | 773
2458906.08654 | 1.27 | 1.38 | 790
2458906.09648 | -2.88 | 1.45 | 837
2458906.10657 | -5.39 | 1.44 | 818

Note. — A machine readable version is available.

##  2 Observations

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 ($\sim 2.7$ hr) with a
$\sim 1$ hour baseline on either side. We used the “C2” decker
($14^{\prime\prime}\times 0\farcs 861$, $R=45,000$) and integrated until the
exposure meter reached 60,000 counts (signal-to-noise ratio (S/N) $\sim 100$
per reduced pixel, $\lesssim 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).

Table 2: Adopted parameters of the WASP-107 System Parameter | Value | Unit | Source
---|---|---|---
$P_{b}$ | $5.7214742$ | days | 1
$t_{c}$ | $7584.329897\pm 0.000032$ | JDaaDays since JD 2,450,000. Sources: (1) Dai & Winn (2017); (2) Piaulet et al. (2021). | 1
$b$ | $0.07\pm 0.07$ |  | 1
$i_{\text{orb},b}$ | $89.887^{+0.074}_{-0.097}$ | degrees | 1
$R_{p}/R_{\star}$ | $0.14434\pm 0.00018$ |  | 1
$a/R_{\star}$ | $18.164\pm 0.037$ |  | 1
$e_{b}$ | $0.06\pm 0.04$ |  | 2
$\omega_{b}$ | $40^{+40}_{-60}$ | degrees | 2
$M_{b}$ | $30.5\pm 1.7$ | M⊕ | 2
$P_{c}$ | $1088^{+15}_{-16}$ | days | 2
$e_{c}$ | $0.28\pm 0.07$ |  | 2
$\omega_{c}$ | $-120^{+30}_{-20}$ | degrees | 2
$M_{c}\sin i_{\text{orb},c}$ | $0.36\pm 0.04$ | MJ | 2
$T_{\text{eff}}$ | $4245\pm 70$ | K | 2
$M_{\ast}$ | $0.683^{+0.017}_{-0.016}$ | M⊙ | 2
$R_{\ast}$ | $0.67\pm 0.02$ | R⊙ | 2
$u_{1}$ | $0.6666\pm 0.0062$ |  | 1
$u_{2}$ | $0.0150\pm 0.0110$ |  | 1

Visually inspecting the observations (Fig. 1) shows an anomalous blueshift
following the transit ingress, followed by an anomalous redshift after the
transit midpoint,111Propagating the uncertainty in $t_{c}$ in Table 2 the
transit midpoint on the night of observation is uncertain to about 9 s.,
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_{\star}$. The expected RM amplitude is $v\sin
i_{\star}(R_{p}/R_{\star})^{2}\sim 40$ m s-1, using previous estimates of
$R_{p}/R_{\star}=0.144$ (Dai & Winn, 2017) and $v\sin i_{\star}\sim 2$ km
$s^{-1}$ (e.g., Anderson et al., 2017). The signal we detected with HIRES is
only $\sim 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_{\star}$ than was spectroscopically inferred (see
Section 3.5), a near-polar orbit, or both.

##  3 Analysis

###  3.1 Rossiter–McLaughlin Model

We used a Gaussian likelihood for the RV time series $(\bm{t},\,\bm{v}_{r})$
given the model parameters $\bm{\Theta}$, and included a RV jitter term
($\sigma_{j}$) to account for additional astrophysical or instrumental noise,

$p(\bm{v}_{r},\,\bm{t}|\bm{\Theta})=\prod_{i=1}^{N}\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left[-\frac{(v_{r,i}-f(t_{i},\,\bm{\Theta}))^{2}}{2\sigma_{i}^{2}}\right],$
(1)

where $\sigma_{i}^{2}=\sigma_{\text{RV,i}}^{2}+\sigma_{j}^{2}$. The model
$f(t_{i},\,\bm{\Theta})$ is given by

$f(t_{i},\,\bm{\Theta})=\text{RM}(t_{i},\,\bm{\theta})+\gamma+\dot{\gamma}(t_{i}-t_{0}),$
(2)

where $\bm{\Theta}=(\bm{\theta},\,\gamma,\,\dot{\gamma})$ is the RM model
parameters ($\bm{\theta}$) as well as an offset ($\gamma$) 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_{0}$ is the time of the first
observation (BJD).

$\text{RM}(t_{i},\,\bm{\theta})$ 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-1, smaller than the jitter from the
Keplerian fit to the full RV sample of $3.9_{-0.4}^{+0.5}$ 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 sky-projected obliquity
($\lambda$), stellar inclination angle ($i_{\star}$), and projected rotational
velocity ($v\sin i_{\star}$). To first order, the impact parameter $b$ and
sky-projected obliquity $\lambda$ determine the shape of the RM signal, while
$v\sin i_{\star}$ and $R_{p}/R_{\star}$ 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 $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}$, $\log(|v_{cb}|)$, $\gamma$,
$\dot{\gamma}$, and $\sigma_{j}$. We placed a uniform hard-bounded prior on
$v\sin i_{\star}\in[0,\,5]$ km s-1 and on $\cos i_{\star}\in[0,\,1]$, and used
a Jeffrey’s prior for $\sigma_{j}$. All other parameters were assigned uniform
priors.

###  3.2 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
$\gamma_{\text{lw}}$, the intrinsic line width, $\zeta$, the line width due to
macroturbulence, given by the Valenti & Fischer (2005) scaling relation

$\zeta=\left(3.98+\frac{T_{\text{eff}}-5770~{}\text{K}}{650~{}\text{K}}\right)~{}\text{km
s}^{-1},$ (3)

and $\beta$, given by

$\beta=\sqrt{\frac{2k_{B}T_{\text{eff}}}{\mu}+\xi^{2}+\beta_{\text{IP}}},$ (4)

where $\xi$ is the dispersion due to microturbulence and $\beta_{\text{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
$\gamma_{\text{lw}},\,\xi$, and $\zeta$ 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 $\lambda$ or $v\sin i_{\star}$. Because of this, we opted to instead adopt
fixed nominal values of $\xi=0.7$ km s-1, $\gamma_{\text{lw}}=1$ km s-1, and
$\zeta=1.63$ km s-1 (from Eq. 3 using $T_{\text{eff}}$ from Table 2).

###  3.3 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 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 $\sim 2$ m s-1, which is significant
given the small amplitude of the RM signal we observe for WASP-107b ($\sim
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
$\lambda$ and $v\sin i_{\star}$ posteriors, we are able to rule out
$|v_{cb}|>450$ m s-1 at 99% confidence, and $>250$ m s-1 at 95% confidence.

Figure 1: The RM effect for WASP-107b. The dark shaded bands show the
16th–84th (black) and 5th–95th (gray) percentiles from the posterior
distribution of the modeled RV. The red best-fit line is the maximum
a-posteriori (MAP) model. The three vertical dashed lines denote, in
chronological order, the times of transit ingress, midpoint, and egress. The
residuals show the data minus the best-fit model. Data points are drawn with
the measurement errors and the best-fit jitter added in quadrature. Figure 2:
Posterior distribution for $\lambda$ and $v\sin i_{\star}$. Although a more
anti-aligned configuration is consistent with the data if $v\sin i_{\star}$ is
small, the most likely orientations are close to polar. A prograde orbit
($|\lambda|<90^{\circ}$) is strongly ruled out. Figure 3: Sky-projected
orbital configuration of WASP-107b’s orbit relative to the stellar rotation
axis. The black lines correspond to posterior draws while the red line is the
MAP orbit from Fig. 1. The direction of WASP-107b’s orbit is denoted by the
red arrow. The stellar rotation axis (black arrow) and lines of stellar
latitude and longitude are drawn for an inclination of $i_{\star}=25^{\circ}$.
The posterior for $i_{\star}$ is illustrated by the shaded gray strip with a
transparency proportional to the probability.

###  3.4 Evidence for a Retrograde/Polar Orbit

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 | aa $v\sin i_{\star}$ is in km s-1 and $v_{cb}$ is in m s-1.
$\sqrt{v\sin i_{\star}}\sin\lambda$ | $-0.126_{-0.771}^{+0.808}$ | -0.72 | aa $v\sin i_{\star}$ is in km s-1 and $v_{cb}$ is in m s-1.
$\cos i_{s}$ | $-0.003_{-0.681}^{+0.682}$ | -0.56 |
$\gamma$ | $0.80_{-1.38}^{+1.36}$ | 0.97 | m s-1
$\dot{\gamma}$ | $-20.83_{-10.94}^{+11.05}$ | -21.85 | m s-2
$\sigma_{\text{jit}}$ | $2.61_{-0.51}^{+0.64}$ | 2.20 | m s-1
$\log(|v_{cb}|)$ | $0.89_{-1.27}^{+1.18}$ | 2.17 | aa $v\sin i_{\star}$ is in km s-1 and $v_{cb}$ is in m s-1.
Derived Parameters
$|\lambda|$ | $118.1_{-19.1}^{+37.8}$ | 112.63 | degrees
$v\sin i_{\star}$ | $0.45_{-0.23}^{+0.72}$ | 0.61 | km s-1
$v_{cb}$ | $-7.74_{-109.71}^{+7.33}$ | -149.41 | m s-1
$i_{\star}$ | $28.17_{-20.04}^{+40.38}$ | 7.06 | degrees
$|\psi|$ | $109.81_{-13.64}^{+28.17}$ | 92.60 | degrees

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 $\lambda$
and $v\sin i_{\star}$ are shown in Fig. 2. A prograde ($|\lambda|<90^{\circ}$)
orbit is ruled out at $>99\%$ confidence. An anti-aligned
($135^{\circ}<\lambda<225^{\circ}$) orbit is allowed if $v\sin i_{\star}$ is
small ($0.26\pm 0.10$ km s-1), although a more polar aligned (but still
retrograde) orbit with $90^{\circ}<|\lambda|<135^{\circ}$ is more likely (if
$v\sin i_{\star}\in[0.22,\,2.09]$ km s-1, 90% CI). The true obliquity $\psi$
will always be closer to a polar orientation than $\lambda$, since $\lambda$
represents the minimum obliquity in the case where the star is viewed edge-on
($i_{\star}=90^{\circ}$). While an equatorial orbit that transits requires
$i_{\star}\sim 90^{\circ}$, 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 $\lambda$ 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 ($\sim 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 $\sim 10\%$ changes on a $\sim 2\%$ transit depth, amounting
to an overall spot depth of $\sim 0.2\%$. Given our estimate of $v\sin
i_{\star}\sim 0.5$ km s-1 this suggests a spot-crossing event would produce a
$\sim 1$ m s-1 RV anomaly, small compared to our measurement uncertainties
($\sim 1.5$ m s-1) and the estimated stellar jitter ($\sim 2.6$ m s-1). In
other words, there is a roughly 33% chance that a spot-crossing event
introduced an additional $0.5\sigma$ 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.

###  3.5 Constraints on the Stellar Inclination

Figure 4: Obliquity of WASP-107b. The true obliquity $\psi$ is calculated
using the constraints on the stellar inclination as inferred from the $v\sin
i_{\star}$ posterior (Section 3.5).

Given a constraint on $v\sin i_{\star}$ and $v$, we can constrain the stellar
inclination $i_{\star}$. 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\pm 0.8$ km s-1, whereas John Brewer (private communication) obtained a
value of $1.5\pm 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_{\star}$ of $<2$ km s-1, as this technique is limited by the
HIRES PSF. 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-1, 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\pm 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\pm 0.02~{}R_{\odot}$ inferred
from the HIRES spectrum (Piaulet et al., 2021) using Specmatch-Emp (Yee et
al., 2017) to constrain the tangential rotational velocity $v=2\pi
R_{\star}/P_{\text{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\leq v\sin i_{\star}$. We
found that $i_{\star}=25.8_{-15.4}^{+22.5}$ degrees (MAP value $7.1^{\circ}$),
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 $\lambda$ near $180^{\circ}$ (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^{\circ}-i_{\star}$) compared to that of our Sun, which has nearly all of
its sunspots contained within $\pm 30^{\circ}$ latitude.

Knowledge of the stellar inclination $i_{\star}$, the orbital inclination
$i_{\text{orb}}$, and the sky-projected obliquity $\lambda$ allows one to
compute the true obliquity $\psi$, as these four angles are related by

$\cos\psi=\cos i_{\text{orb}}\cos i_{\star}+\sin i_{\text{orb}}\sin
i_{\star}\cos\lambda.$ (5)

The resulting posterior distribution for the true obliquity $\psi$ 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 $\lambda$, due to the nearly pole-on
viewing geometry of the star itself.

##  4 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}\sqrt{a_{b}}\ll M_{c}\sqrt{a_{c}}$), and since
$a_{b}/a_{c}\ll 1$, we can approximate the Hamiltonian by expanding to
quadrupole order in semimajor axis ratio

$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\frac{1}{16}n_{b}\frac{M_{c}}{M_{\star}}\left(\frac{a_{b}}{a_{c}\sqrt{1-e_{c}^{2}}}\right)^{3}\left[\frac{(5-3G_{b}^{2})(3H_{b}^{2}-G_{b}^{2})}{G_{b}^{2}}\right.$
(6)
$\displaystyle\left.+\frac{15(1-G_{b}^{2})(G_{b}^{2}-H_{b}^{2})\cos(2g_{b})}{G_{b}^{2}}\right]+\frac{GM_{\star}}{a_{b}c^{2}}\frac{3n_{b}}{G_{b}},$

where the last term is the addition from general relativity (GR) and
$n_{b}=2\pi/P_{b}$. The quantities $G$ and $H$ are the canonical Delaunay
variables

$\displaystyle G_{b}$ $\displaystyle=\sqrt{1-e_{b}^{2}}$
$\displaystyle\leftrightarrow\;g_{b}=\omega_{b},$ (7) $\displaystyle H_{b}$
$\displaystyle=G\cos i_{b}$ $\displaystyle\leftrightarrow\;h_{b}=\Omega_{b},$

where the double-arrow ($\leftrightarrow$) symbolizes conjugate variables,
$\omega_{b}$ is the argument of perihelion of the inner planet, $\Omega_{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 $\propto Ma^{1/2}$). With this approximation, $i_{b}$ is
the relative inclination between the two planets.

Figure 5: Evolution of WASP-107b’s true obliquity ($\psi_{b}$, solid line)
throughout the the $N$-body simulation using the system parameters given in
Table 2. The outer planet has $M_{c}=M\sin i_{\text{orb},c}$ and was
initialized with an obliquity of $\psi_{c}=60^{\circ}$ (dashed line). The
obliquity of planet b oscillates between $\psi_{c}\pm\psi_{c}$ every $\sim
2.5$ Myr due to nodal precession. If $\sin i_{\text{orb},c}<1$ then the larger
$M_{c}$ simply produces a shorter nodal precession timescale. The right panel
shows the evolution of the inclinations with the difference in the longitudes
of ascending node.

###  4.1 Kozai–Lidov oscillations

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 $\sim 39\fdg 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 $\pi$ Men (Xuan & Wyatt, 2020; Yee et al., 2018).
The precession rate from GR is given by

$\dot{\omega}_{GR}=\frac{GM_{\star}}{a_{b}c^{2}}\frac{3n_{b}}{G_{b}^{2}},$ (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

$\tau_{\text{Kozai}}=\frac{2P_{c}^{2}}{3\pi
P_{b}^{2}}\frac{M_{\star}}{M_{c}}(1-e_{c}^{2})^{3/2}\approx
210,000~{}\text{yr},$ (9)

five times longer. The condition for Kozai–Lidov cycles to be suppressed by
relativistic precession is $\tau_{\text{Kozai}}\dot{\omega}_{\text{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
$\tau_{\text{Kozai}}$ and $\tau_{\text{GR}}$. While the true mass of WASP-107c
is likely to be larger than the derived $M\sin i_{\text{orb},c}$, it would
need to be $\sim 10$ times larger for Kozai–Lidov oscillations to occur. This
would imply a near face-on orbit of at most $i_{\text{orb},c}<5\fdg 5$. 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.

###  4.2 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 $\pi$ 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 $\psi_{\text{c}}\sim 55^{\circ}$. Then the longitude of ascending node
$\Omega_{b}$ evolves in a secular manner according to Yee et al. (2018),

$\frac{d\Omega_{b}}{dt}=\frac{\partial\mathcal{H}}{\partial
H_{b}}=\frac{n_{b}}{8}\frac{M_{c}}{M_{\star}}\left(\frac{a_{b}}{a_{c}\sqrt{1-e_{c}^{2}}}\right)^{3}\left(\frac{15-9G_{b}^{2}}{G_{b}^{2}}\right)H_{b}.$
(10)

The associated timescale $\tau_{\Omega_{b}}=2\pi/\dot{\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 $\approx\psi_{c}\pm\psi_{c}$. Thus at certain
times the observer may see a highly misaligned orbit ($\psi_{b}\sim
2\psi_{c}$) for the inner planet, while at other times the observer may see an
aligned orbit ($\psi_{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, ${\sim}2\psi_{c}=120^{\circ}$) 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 $\psi_{b}$ oscillates in the range
$0^{\circ}\text{--}120^{\circ}$ due to the precession of $\Omega_{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 $\gtrsim 80\%$
the amplitude from nodal precession (${\sim}2\psi_{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 $\psi_{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 $\pi$ 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 $\pi$ Men systems (Xuan & Wyatt, 2020), although the inner
planet in $\pi$ Men is only slightly misaligned with $\lambda=24\pm 4.1$
degrees (Hodžić et al., submitted), 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 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.

Figure 6: Top: polar plot showing the absolute sky-projected obliquity as the
azimuthal coordinate and normalized orbital distance as the radial coordinate,
for ${<}100~{}M_{\oplus}$ planets around stars with $T_{\text{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 $\pi$ 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 $\psi$ is assumed to vary as
$\cos[(\pi/2)\psi(t)/\psi_{\text{max}}]=\sin^{2}(\pi t/\tau)$, where
$t\in[0,\tau=1]$. This recreates the shape of the oscillating inclination in
Fig. 5. The amplitude $\psi_{\text{max}}$ is twice the outer planet’s
inclination which is plotted for three different distributions (shown on the
right): uniform between $[0^{\circ},\,90^{\circ}]$ (gray), uniform between
$[40^{\circ},\,60^{\circ}]$ (red), and using the von-Mises Fisher distribution
from Masuda et al. (2020) calculated in a hierarchical manner incorporating
their posterior distribution for the shape parameter $\sigma$ 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_{\text{orb},b}=90^{\circ}$) and is marginalized over stellar inclination
angle (solid histogram). $M_{p}<100~{}M_{\oplus}$ 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 $\sim 40$–$60^{\circ}$.

###  4.3 Disk dispersal-driven tilting

Recently, Petrovich et al. (2020) showed that, even for $\psi_{c}\sim
0^{\circ}$, 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\mapsto c$ and $c\mapsto\text{disk}$). The outer
planet still induces a nodal precession on the inner planet according to Eq.
10. If at first the rate $d\Omega_{c}/dt>d\Omega_{b}/dt$, then as the disk
dissipates (and $M_{\text{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^{\circ}$. 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 ($\eta_{\text{GR}}$) with the stellar quadrupole moment
($\eta_{\star}$). We found that $\eta_{\text{GR}}>\eta_{\star}+6$ at 99.76%
confidence, $\eta_{\star}+6>\eta_{\text{GR}}>4$ at 0.155% confidence, and
$\eta_{\text{GR}}<4$ at 0.084% confidence (i.e., $\eta_{\text{GR}}\sim 30-80$
and $\eta_{\star}\sim 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_{\text{orb},c}$ and $P_{c}$ and integrating over all possible initial
obliquities for the outer planet. Evaluating their Eq. (3), we found that the
resonance that drives the inner planet to high obliquities is always crossed.
We calculated the adiabatic parameter
$x_{\text{ad}}\equiv\tau_{\text{disk}}/\tau_{\text{adia}}$ from the disk
dispersal timescale and the adiabatic time (their Eq. 7), taking
$\tau_{\text{disk}}$ to be 1 Myr. In the orbital configurations where
$x_{\text{ad}}>1$ (adiabatic crossing) we computed the final obliquity from
their Eq. (12) ($I_{\text{crit}}$). Otherwise, the final obliquity was set to
$I_{\text{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 $\sim 8^{\circ}$, and
is guaranteed for $\psi_{\text{init},c}\gtrsim 25^{\circ}$. 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~{}\text{AU})^{3}>(M_{c}/0.5~{}M_{J})$. Since we only have a
constraint on $M\sin i_{\text{orb},c}$, this condition is satisfied if
$i_{\text{orb},c}\in[60^{\circ}-90^{\circ}]$. Such a viewing geometry, in
conjunction with an obliquity of $\psi_{c}>25^{\circ}$, 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 $\psi=90^{\circ}$ polar
orbit from below and stops at $\psi_{b}=90^{\circ}$. 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 ($\psi_{c}\gtrsim
40-50^{\circ}$) obliquity needed to excite either Kozai–Lidov or nodal
precession cycles.

##  5 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_{\star}$. This low $v\sin i_{\star}$ implies that we are
viewing the star close to one of its poles, reinforcing the near-polar 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 orbital 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, in preparation) are consistent with a circular orbit. The eccentricity
damping timescale due to tidal forces is only $\sim 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 $\pi$ 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 $\mu$as), prevents
a detection of the outer planet using astrometry. Assuming future Gaia data
releases have the same astrometric precision as in DR2 ($44~{}\mu$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. 2020; and 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.

We thank Konstantin Batygin, Cristobol Petrovich, and Jerry Xuan for helpful
comments and productive discussions on orbital dynamics, and Josh Winn for
constructive feedback that improved this manuscript. R.A.R. and A.C.
acknowledge support from the National Science Foundation through the Graduate
Research Fellowship Program (DGE 1745301, DGE 1842402). C.D.D. acknowledges
the support of the Hellman Family Faculty Fund, the Alfred P. Sloan
Foundation, the David & Lucile Packard Foundation, and the National
Aeronautics and Space Administration via the TESS Guest Investigator Program
(80NSSC18K1583). I.J.M.C. acknowledges support from the NSF through grant
AST-1824644. D.H. acknowledges support from the Alfred P. Sloan Foundation,
the National Aeronautics and Space Administration (80NSSC18K1585,
80NSSC19K0379), and the National Science Foundation (AST-1717000). E.A.P.
acknowledges the support of the Alfred P. Sloan Foundation. L.M.W. is
supported by the Beatrice Watson Parrent Fellowship and NASA ADAP Grant
80NSSC19K0597. We thank the time assignment committees of the University of
California, the California Institute of Technology, NASA, and the University
of Hawai‘i for supporting the TESS–Keck Survey with observing time at the W.
M. Keck Observatory. We gratefully acknowledge the efforts and dedication of
the Keck Observatory staff for support of HIRES and remote observing. We
recognize and acknowledge the cultural role and reverence that the summit of
Maunakea has within the indigenous Hawaiian community. We are deeply grateful
to have the opportunity to conduct observations from this mountain. Keck I
(HIRES)

## References

  * Allart et al. (2019) Allart, R., Bourrier, V., Lovis, C., et al. 2019, A&A, 623, A58
  * Anderson et al. (2017) Anderson, D. R., Collier Cameron, A., Delrez, L., et al. 2017, A&A, 604, A110
  * Bourrier et al. (2017) Bourrier, V., Cegla, H. M., Lovis, C., & Wyttenbach, A. 2017, A&A, 599, A33
  * Bourrier et al. (2018) Bourrier, V., Lovis, C., Beust, H., et al. 2018, Nature, 553, 477
  * Brewer et al. (2016) Brewer, J. M., Fischer, D. A., Valenti, J. A., & Piskunov, N. 2016, ApJS, 225, 32
  * Butler et al. (1996) Butler, R. P., Marcy, G. W., Williams, E., et al. 1996, PASP, 108, 500
  * Cubillos et al. (2017) Cubillos, P., Harrington, J., Loredo, T. J., et al. 2017, AJ, 153, 3
  * Dai & Winn (2017) Dai, F., & Winn, J. N. 2017, AJ, 153, 205
  * Dai et al. (2020) Dai, F., Roy, A., Fulton, B., et al. 2020, AJ, 160, 193
  * Dalal et al. (2019) Dalal, S., Hébrard, G., Lecavelier des Étangs, A., et al. 2019, A&A, 631, A28
  * Dalba et al. (2020) Dalba, P. A., Gupta, A. F., Rodriguez, J. E., et al. 2020, AJ, 159, 241
  * De Rosa et al. (2020) De Rosa, R. J., Dawson, R., & Nielsen, E. L. 2020, A&A, 640, A73
  * Esposito, M. et al. (2014) Esposito, M., Covino, E., Mancini, L., et al. 2014, A&A, 564, L13
  * Fabrycky & Tremaine (2007) Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298
  * Ford (2006) Ford, E. B. 2006, ApJ, 642, 505
  * Foreman-Mackey (2016) Foreman-Mackey, D. 2016, JOSS, 1
  * Foreman-Mackey et al. (2013) Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306
  * Gelman et al. (2003) Gelman, A., Carlin, J. B., Stern, H. S., & Rubin, D. B. 2003, Bayesian Data Analysis, 2nd edn. (Chapman and Hall)
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357–362
  * Hébrard et al. (2008) Hébrard, G., Bouchy, F., Pont, F., et al. 2008, A&A, 488, 763
  * Hirano et al. (2011) Hirano, T., Suto, Y., Winn, J. N., et al. 2011, ApJ, 742, 69
  * Hodžić et al. (submitted) Hodžić, V. K., Triaud, A. H. M. J., Cegla, H. M., Chaplin, W. J., & Davies, G. R. submitted, MNRAS
  * Høg et al. (2000) Høg, E., Fabricius, C., Makarov, V. V., et al. 2000, A&A, 355, L27
  * Howard et al. (2010) Howard, A. W., Johnson, J. A., Marcy, G. W., et al. 2010, ApJ, 721, 1467
  * Howell et al. (2014) Howell, S. B., Sobeck, C., Haas, M., et al. 2014, PASP, 126, 398
  * Huang et al. (2017) Huang, C. X., Petrovich, C., & Deibert, E. 2017, AJ, 153, 210
  * Huber et al. (2013) Huber, D., Carter, J. A., Barbieri, M., et al. 2013, Science, 342, 331
  * Hunter (2007) Hunter, J. D. 2007, CSE, 9, 90
  * Kamiaka et al. (2019) Kamiaka, S., Benomar, O., Suto, Y., et al. 2019, The Astronomical Journal, 157, 137
  * Kirk et al. (2020) Kirk, J., Alam, M. K., López-Morales, M., & Zeng, L. 2020, AJ, 159, 115
  * Kiseleva et al. (1998) Kiseleva, L. G., Eggleton, P. P., & Mikkola, S. 1998, MNRAS, 300, 292
  * Kostov et al. (2016) Kostov, V. B., Moore, K., Tamayo, D., Jayawardhana, R., & Rinehart, S. A. 2016, ApJ, 832, 183
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591
  * Kreidberg et al. (2018) Kreidberg, L., Line, M. R., Thorngren, D., Morley, C. V., & Stevenson, K. B. 2018, ApJL, 858, L6
  * Lidov (1962) Lidov, M. 1962, Planet. Space Sci., 9, 719
  * Mann et al. (2020) Mann, A. W., Johnson, M. C., Vanderburg, A., et al. 2020, AJ, 160, 179
  * Mardling (2010) Mardling, R. A. 2010, MNRAS, 407, 1048
  * Masuda & Winn (2020) Masuda, K., & Winn, J. N. 2020, AJ, 159, 81
  * Masuda et al. (2020) Masuda, K., Winn, J. N., & Kawahara, H. 2020, AJ, 159, 38
  * McLaughlin (1924) McLaughlin, D. B. 1924, ApJ, 60, 22
  * Millholland et al. (2020) Millholland, S., Petigura, E., & Batygin, K. 2020, ApJ, 897, 7
  * Montet et al. (2020) Montet, B. T., Feinstein, A. D., Luger, R., et al. 2020, AJ, 159, 112
  * Močnik et al. (2017) Močnik, T., Hellier, C., Anderson, D. R., Clark, B. J. M., & Southworth, J. 2017, MNRAS, 469, 1622
  * Oklopčić & Hirata (2018) Oklopčić, A., & Hirata, C. M. 2018, ApJL, 855, L11
  * Palle et al. (2020) Palle, E., Oshagh, M., Casasayas-Barris, N., et al. 2020, A&A, 643, A25
  * Petrovich et al. (2020) Petrovich, C., Muñoz, D. J., Kratter, K. M., & Malhotra, R. 2020, ApJL, 902, L5
  * Piaulet et al. (2021) Piaulet, C., Benneke, B., Rubenzahl, R. A., et al. 2021, AJ, 161, 70
  * Powell (1964) Powell, M. J. D. 1964, The Computer Journal, 7, 155
  * Queloz et al. (2000) Queloz, D., Eggenberger, A., Mayor, M., et al. 2000, A&A, 359, L13
  * Queloz et al. (2010) Queloz, D., Anderson, D., Collier Cameron, A., et al. 2010, A&A, 517, L1
  * Rein & Liu (2012) Rein, H., & Liu, S.-F. 2012, A&A, 537, A128
  * Rein & Tamayo (2015) Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376
  * Rossiter (1924) Rossiter, R. A. 1924, ApJ, 60, 15
  * Sanchis-Ojeda et al. (2013) Sanchis-Ojeda, R., Winn, J. N., Marcy, G. W., et al. 2013, ApJ, 775, 54
  * Shporer & Brown (2011) Shporer, A., & Brown, T. 2011, ApJ, 733, 30
  * Southworth (2011) Southworth, J. 2011, MNRAS, 417, 2166
  * Spake et al. (2018) Spake, J. J., Sing, D. K., Evans, T. M., et al. 2018, Nature, 557, 68
  * Spalding & Batygin (2015) Spalding, C., & Batygin, K. 2015, ApJ, 811, 82
  * Stefánsson et al. (2020) Stefánsson, G., Kopparapu, R., Lin, A., et al. 2020, AJ, 160, 259
  * Tamayo et al. (2019) Tamayo, D., Rein, H., Shi, P., & Hernandez, D. M. 2019, MNRAS, 491, 2885
  * Valenti & Fischer (2005) Valenti, J. A., & Fischer, D. A. 2005, ApJSupplement Series, 159, 141
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, NatMe, 17, 261
  * Vogt et al. (1994) Vogt, S. S., Allen, S. L., Bigelow, B. C., et al. 1994, in Proc. SPIE, Vol. 2198, Instrumentation in Astronomy VIII, ed. D. L. Crawford & E. R. Craine, 362
  * Winn et al. (2010a) Winn, J. N., Fabrycky, D., Albrecht, S., & Johnson, J. A. 2010a, ApJL, 718, L145
  * Winn et al. (2010b) Winn, J. N., Johnson, J. A., Howard, A. W., et al. 2010b, ApJL, 723, L223
  * Xuan & Wyatt (2020) Xuan, J. W., & Wyatt, M. C. 2020, MNRAS, 497, 2096
  * Yee et al. (2017) Yee, S. W., Petigura, E. A., & von Braun, K. 2017, ApJ, 836, 77
  * Yee et al. (2018) Yee, S. W., Petigura, E. A., Fulton, B. J., et al. 2018, AJ, 155, 255
  * Zhou et al. (2020) Zhou, G., Winn, J. N., Newton, E. R., et al. 2020, ApJ, 892, L21