# HATS-38 b and WASP-139 b join a growing group of eccentric hot Neptunes on
polar orbits111Based on observations made with ESO Telescopes at the La Silla
Paranal Observatory under programs ID 111.24VT.001, 111.24VT.002, and
112.25W1.001

Juan I. Espinoza-Retamal Instituto de Astrofísica, Pontificia Universidad
Católica de Chile, Av. Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile
Millennium Institute for Astrophysics, Santiago, Chile  Anton Pannekoek
Institute for Astronomy, University of Amsterdam, Science Park 904, 1098 XH
Amsterdam, The Netherlands  Guðmundur Stefánsson Anton Pannekoek Institute for
Astronomy, University of Amsterdam, Science Park 904, 1098 XH Amsterdam, The
Netherlands  Cristobal Petrovich Instituto de Astrofísica, Pontificia
Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 782-0436 Macul,
Santiago, Chile  Millennium Institute for Astrophysics, Santiago, Chile
Department of Astronomy, Indiana University, Bloomington, IN 47405, USA
Rafael Brahm Millennium Institute for Astrophysics, Santiago, Chile  Facultad
de Ingeniería y Ciencias, Universidad Adolfo Ibáñez, Av. Diagonal las Torres
2640, Peñalolén, Santiago, Chile  Data Observatory Foundation, Santiago, Chile
Andrés Jordán Millennium Institute for Astrophysics, Santiago, Chile  Facultad
de Ingeniería y Ciencias, Universidad Adolfo Ibáñez, Av. Diagonal las Torres
2640, Peñalolén, Santiago, Chile  Data Observatory Foundation, Santiago, Chile
El Sauce Observatory — Obstech, Coquimbo, Chile  Elyar Sedaghati Millennium
Institute for Astrophysics, Santiago, Chile  Facultad de Ingeniería y
Ciencias, Universidad Adolfo Ibáñez, Av. Diagonal las Torres 2640, Peñalolén,
Santiago, Chile  European Southern Observatory (ESO), Av. Alonso de Córdova
3107, 763 0355 Vitacura, Santiago, Chile  Jennifer P. Lucero Instituto de
Astrofísica e Ciências do Espaço, Universidade do Porto, CAUP, Rua das
Estrelas, 4150-762 Porto, Portugal  Departamento de Física e Astronomia,
Faculdade de Ciências, Universidade do Porto, Rua do Campo Alegre, 4169-007
Porto, Portugal  Marcelo Tala Pinto Millennium Institute for Astrophysics,
Santiago, Chile  Facultad de Ingeniería y Ciencias, Universidad Adolfo Ibáñez,
Av. Diagonal las Torres 2640, Peñalolén, Santiago, Chile  Diego J. Muñoz
Department of Astronomy and Planetary Science, Northern Arizona University,
Flagstaff, AZ 86011, USA  Gavin Boyle El Sauce Observatory — Obstech,
Coquimbo, Chile  Cavendish Laboratory, J J Thomson Avenue, Cambridge, CB3 0HE,
UK  Rodrigo Leiva Millennium Institute for Astrophysics, Santiago, Chile
Instituto de astrofísica de Andalucía, CSIC, Glorieta de la Astronomía s/n,
18008 Granada, Spain  Vincent Suc Millennium Institute for Astrophysics,
Santiago, Chile  Facultad de Ingeniería y Ciencias, Universidad Adolfo Ibáñez,
Av. Diagonal las Torres 2640, Peñalolén, Santiago, Chile  El Sauce Observatory
— Obstech, Coquimbo, Chile

###### Abstract

We constrain the sky-projected obliquities of two low-density hot Neptune
planets, HATS-38 b and WASP-139 b, orbiting nearby G and K stars using
Rossiter-McLaughlin (RM) observations with VLT/ESPRESSO, yielding
$\lambda=-108_{-16}^{+11}$ deg and $-85.6_{-4.2}^{+7.7}$ deg, respectively. To
model the RM effect, we use a new publicly available code, ironman, which is
capable of jointly fitting transit photometry, Keplerian radial velocities,
and RM effects. The two planets have residual eccentricities
($e=0.112_{-0.070}^{+0.072}$, and $0.103_{-0.041}^{+0.050}$, respectively),
and together with the obliquity constraints, we show that they join a growing
group of eccentric hot and low-density Neptunes on polar orbits. We use long-
term radial velocities to rule out companions with masses $\sim 0.3-50$
$M_{J}$ within $\sim$10 au. We show that the orbital architectures of the two
Neptunes disfavor an origin from primordial disk misalignment and/or in-situ
secular interactions with distant companions and instead favor high-
eccentricity migration from $\gtrsim 2$ au driven by a distant companion.
Finally, we performed a hierarchical Bayesian modeling of the true obliquity
distribution of Neptunes and found suggestive evidence for a higher
preponderance of polar orbits of hot Neptunes compared to Jupiters. However,
we note that the exact distribution is sensitive to the choice of priors,
highlighting the need for additional obliquity measurements of Neptunes to
robustly compare the hot Neptune obliquity distribution to Jupiters.

Exoplanets (498) — Hot Neptunes (754) — Exoplanet dynamics (490) — Planetary
alignment (1243)

##  1 Introduction

More than 5,000 exoplanets have been confirmed, with super-Earths and mini-
Neptunes being the most prevalent. Essential information about the formation
of these planetary systems can be found in their architectures, as they serve
as a powerful probe of their dynamical history. One of the main signatures of
the dynamical evolution of systems is the stellar obliquity ($\psi$)—the angle
between the stellar spin axis and the planet’s orbital axis. Significant
progress has been made in recent years in measuring the sky-projected
obliquity ($\lambda$) distribution of hot Jupiters. Measurements performed via
the Rossiter-McLaughlin (RM) effect have uncovered a broad distribution of
$\lambda$ for hot Jupiters, ranging from well-aligned systems to systems on
highly misaligned orbits (Albrecht et al., 2012). These results have been
interpreted as evidence that hot Jupiter formation involves dynamical
perturbations that excite orbital inclinations.

However, the signature of the dynamical evolution is the true 3D obliquity
$\psi$, not its sky projection $\lambda$. The true obliquity $\psi$, can be
estimated from $\lambda$ measured by the RM effect through combining it with
constraints on the orbital inclination, $i$ constrained by the shape of a
planetary transit, and the stellar inclination, $i_{\star}$, which can be
constrained by measuring the stellar rotation period ($P_{\rm rot}$), the
projected rotational velocity ($v\sin i_{\star}$), and the radius of the host
star (e.g., Hirano et al., 2014; Morton & Winn, 2014; Masuda & Winn, 2020).

By systematically calculating true obliquities from sky-projected obliquities,
Albrecht et al. (2021) revealed an architectural dichotomy: hot Jupiter orbits
appear to separate into well-aligned ($\psi\sim 0$ deg) and polar ($\psi\sim
90$ deg) orbits. However, more recent work by both Siegel et al. (2023) and
Dong & Foreman-Mackey (2023), using Bayesian approaches, did not find strong
evidence for a dichotomy. In turn, the Bayesian modeling suggests the $\psi$
distribution is unimodal and peaked at $0$ deg with an almost isotropic tail
for misaligned systems, with no significant clustering at $90$ deg. As hot
Jupiters represent an infrequent outcome of planet formation (Batalha et al.,
2013), it is unclear if the observed obliquity distribution is related to the
formation of hot Jupiters, or instead reflects more general aspects of star
and planet formation.

Thanks to recent observational efforts, the current sample of planets with
measured obliquity shows a growing number of smaller planets including hot
Neptunes ($10<M_{p}/M_{\oplus}<50$ or $2<R_{p}/R_{\oplus}<6$) in the ‘hot
Neptune desert’, a region of the period-radius diagram where Neptune-sized
planets are notably rare (Mazeh et al., 2016) likely due to intense stellar
radiation causing atmospheric erosion. Interestingly, many hot Neptunes are
observed to be in close-to-polar orbits and have evaporating atmospheres,
including HAT-P-11 b (Sanchis-Ojeda & Winn, 2011), GJ 436 b (Bourrier et al.,
2018, 2022), WASP-107 b (Dai & Winn, 2017; Rubenzahl et al., 2021), and GJ
3470 b (Stefansson et al., 2022), suggesting a possible link between the
Neptune desert and the dynamical history of low-mass gas giants. Theoretical
scenarios that are capable of explaining the polar orbits of these planets
require the presence of massive distant companions (e.g., Petrovich et al.,
2020). Recently, outer companions have been confirmed in the HAT-P-11 (Yee et
al., 2018) and WASP-107 (Piaulet et al., 2021) systems, suggesting that the
transiting planets arrived at their current polar orbits through dynamical
interactions (see Lu et al., 2024; Yu & Dai, 2024). As for the atmospheres,
the orbital eccentricities may be large enough so that tidal heating can
inflate these planets (e.g., by a factor of $\simeq 1.5-1.7$ for GJ 3470b;
Stefansson et al., 2022), possibly accounting for their evaporation.
Nevertheless, the sample of hot Neptunes with measured obliquities is still
small, and more measurements are necessary to test these emerging trends and
asses the underlying mechanism behind the formation of this transiting
population.

In this work, we present ESPRESSO (Pepe et al., 2021) observations of the RM
effect of the two hot Neptunes HATS-38 b (Jordán et al., 2020) and WASP-139 b
(Hellier et al., 2017). HATS-38 b is a hot Neptune in the middle of the hot
Neptune desert orbiting a G-dwarf host star, while WASP-139 b is a super-
Neptune orbiting a K-star. The precise in-transit spectroscopic observations
with ESPRESSO reveal that both planets are in nearly-polar and eccentric
orbits, continuing with the observed trends for low-density gas giants. We
also performed a long-term radial velocity (RV) monitoring of both stars in
order to constrain the presence of distant planetary companions in the
systems.

We describe our observations in Section 2. In Section 3, we present our
stellar analysis for both stars. In Section 4, we update the orbital ephemeris
of HATS-38 b and WASP-139 b, while in Section 5, we describe the joint
analysis of photometry and RV data to get the parameters of the planets and
their orbits. We discuss the implications of these measurements in Section 6
and summarize our findings in Section 7.

Figure 1: Observations of HATS-38 along with the best-fit model. All error
bars include a white noise jitter term. a) ESPRESSO observations of the RM
effect. The red line and shaded area correspond to the best and $1\sigma$
models, respectively. b) Phase-folded out-of-transit RVs showing $1\sigma$,
$2\sigma$, and $3\sigma$ models. c) Residuals of the model for the RVs as a
function of time. d-j) Different light curves. We show the binned TESS Year 5
data in black.

Figure 2: Same as Figure 1 but for WASP-139. White points in the RM effect
plot were not considered for the fit given their low S/N due to increased
cloud coverage during those exposures. Including or removing those points
results in the same constraints on the sky-projected obliquity. Table 1:
Summary of observations used in this work.

Facility | Date | Notes | Reference
---|---|---|---
HATS-38 |  |  |
ESPRESSO | 2023 Apr 18 | 26 RVs | This work
HARPSa | 2016 Nov–2017 May | 24 RVs | Jordán et al. (2020)
FEROS | 2016 Dec–2017 Mar | 10 RVs | Jordán et al. (2020)
PFS | 2017 Apr 05-19 | 5 RVs | Jordán et al. (2020)
TESSb | Years 1, 3 & 5 | 30, 10 & 2 min cadence | This work
MOANA/ES | 2023 Apr 18 |  $r^{\prime}$ filter | This work
CHAT | 2017 Feb 05 |  $r$ filter | Jordán et al. (2020)
LCO/SAAO | 2017 Mar 30 |  $i^{\prime}$ filter | Jordán et al. (2020)
LCO/CTIO | 2017 Apr 03 |  $i^{\prime}$ filter | Jordán et al. (2020)
WASP-139 |  |  |
ESPRESSO | 2023 August 08 | 22 RVs | This work
HARPSa | 2014 Sep-2015 Jan | 23 RVs | Hellier et al. (2017)
CORALIEc | 2008 Oct–2015 Dec | 24 RVs | Hellier et al. (2017)
TESSb | Years 1 & 3 | 2 min cadence | This work
TRAPPIST | 2014 Aug 06 |  $I+z$ band | Hellier et al. (2017)
TRAPPIST | 2015 Sep 07 |  $I+z$ band | Hellier et al. (2017)
EulerCAM | 2015 Sep 07 |  $NGTS$ filter | Hellier et al. (2017)

Note. — a We took two extra HARPS measurements per target and reduced all data
again in a uniform way (see Section 2.2).b We combined all the TESS sectors of
each year, and in the analysis, we assumed that different years were different
instruments. c We considered the pre and post CORALIE update data as taken by
different instruments.

##  2 Observations

Table 1 shows a summary of the observations used in this work. Next, we
describe novel observations of HATS-38 and WASP-139.

###  2.1 ESPRESSO Transit Spectroscopy

We observed a single transit of each planet, HATS-38 b and WASP-139 b, with
the ESPRESSO spectrograph (Pepe et al., 2021). ESPRESSO is a highly-
stabilized, fiber-fed cross-dispersed echelle spectrograph installed at the
Incoherent Combined Coudé Focus of ESO’s Paranal Observatory in Chile. It
covers a wavelength range from 380 to 788 nm at a resolving power of
$R\approx$ 140,000 in single Unit Telescope (UT) high-resolution mode. The
transit of HATS-38 b was observed on the night of 2023 April 18, between 01:25
and 06:01 UTC. We obtained 26 spectra of the host star during the primary
transit with UT1 at an exposure time of 610 s. The observations were performed
under clear sky conditions, with atmospheric seeing in the range of
$0.42-1.7^{\prime\prime}$. The spectra have a median signal-to-noise ratio
(S/N) of 39 at 550 nm and a median RV uncertainty of 1.7 m/s. The transit of
WASP-139 b was observed on the night of 2023 August 8, between 05:35 and 10:04
UTC. We obtained 22 spectra of the host star during the primary transit with
UT3. The sky was cloudy at the beginning of observations before the transit,
and conditions improved during the night. We started observations with an
exposure time of 900 s, and then changed it to 600 s when conditions improved.
In our analysis, we excluded 4 of the 22 spectra given their low S/N ($<20$ at
550 nm). The remaining 18 spectra have a median S/N of 31 at 550 nm and a
median RV uncertainty of 2.2 m/s. Both datasets were reduced with the ESO
Reflex environment (Freudling et al., 2013) using the dedicated ESPRESSO data
reduction pipeline (v. 3.1.0, Sosnowska et al., 2015; Modigliani et al.,
2020), including all standard reduction steps, which also provides RVs by
fitting a Gaussian model to the calculated cross-correlation function (CCF).
The CCF was calculated at steps of 0.5 km/s (representing the sampling of the
spectrograph) for $\pm 20$ km/s centered on the estimated systemic velocity.
The RVs showing the RM effect of both targets are presented in panels a of
Figures 1 and 2.

###  2.2 HARPS Follow-up Spectroscopy

HATS-38 and WASP-139 were both observed with the HARPS spectrograph mounted on
the 3.6m telescope installed at the ESO La Silla Observatory. HARPS is a high-
resolution ($R\approx$ 120,000) stabilized instrument that uses a secondary
fiber to trace instrumental velocity variations by measuring drifts of a
Fabry-Perot interferometry pattern. Archival observations were available
through the ESO archive. These observations were used initially to confirm the
planetary nature of both planets. Besides these archival HARPS data, we
obtained two additional spectra for each system in order to look for long-
period companions. These new observations were performed in November and
December of 2023. For HATS-38 we used an exposure time of 1500 s while for
WASP-139 we adopted an exposure time of 1800 s. Archival and new HARPS data
were homogeneously processed with the ceres (Brahm et al., 2017a) pipeline,
which performs all processing steps to obtain precision RVs starting from the
raw images.

###  2.3 Observatoire Moana Photometry

Simultaneously with ESPRESSO observations, we observed the transit of HATS-38
b using the station of the Observatoire Moana located in El Sauce (ES)
Observatory in Chile. Observatoire Moana is a global network of small-aperture
robotic optical telescopes. The ES station consists of a 0.6 m CDK robotic
telescope coupled to an Andor iKon-L 936 deep depletion 2k $\times$ 2k CCD
with a scale of 0.67′′ per pixel. For these observations, we used a Sloan
$r^{\prime}$ filter, and the exposure time was set to 80 s. The MOANA/ES light
curve is shown in panel h of Figure 1, along with the best model.

###  2.4 TESS Photometry

For this work, we also used all the available TESS (Ricker et al., 2015) light
curves of both targets. HATS-38 was observed in Sector 9 (Year 1) at a cadence
of 30 minutes, in Sectors 35 and 36 (Year 3) at a cadence of 10 minutes, and
in Sector 62 (Year 5) at a cadence of 2 minutes. WASP-139 was observed by TESS
in Sectors 3 and 4 (Year 1), and 31 and 32 (Year 3) at a cadence of 2 minutes.
We searched and downloaded all the light curves using lightkurve (Lightkurve
Collaboration et al., 2018). For this work, we used TESS light curves that
were processed by the Science Processing Operations Center (SPOC) pipeline
(Jenkins et al., 2016). In particular, we used the Presearch Data Conditioning
(PDC) light curves of the SPOC pipeline, which are corrected for pointing or
focus-related instrumental signatures, discontinuities resulting from
radiation events in the CCD detectors, outliers, and flux contamination. TESS
PDC light curves, together with the best model, are shown in the bottom panels
of Figures 1 and 2.

##  3 Stellar Parameters

We obtained the stellar parameters of HATS-38 and WASP-139 following the
procedure presented in Brahm et al. (2019). In brief, we used a two-step
iterative process.

The first process consists of obtaining the atmospheric parameters ($T_{\rm
eff}$, $\log{g}$, [Fe/H], and $v\sin{i}$) of the host star from a high-
resolution spectrum using the zaspe code (Brahm et al., 2017b). This code
compares the observed spectrum to a grid of synthetic spectra and determines
reliable uncertainties that take into account systematic mismatches between
the observations and the imperfect theoretical models. For this analysis, we
used the co-added out-of-transit HARPS spectra to obtain the atmospheric
parameters.

The second step consists of obtaining the physical parameters of the star by
using publicly available broad-band magnitudes of the star and comparing them
with those produced by different PARSEC stellar evolutionary models (Bressan
et al., 2012) by taking into account the distance to the star computed from
the Gaia DR2 (Gaia Collaboration et al., 2018) parallax. This procedure
delivers a new value of $\log{g}$ that is held fixed in a new run of zaspe. We
iterate between the two procedures until reaching convergence, which happens
when two consecutive zaspe runs deliver the same values of $T_{\rm eff}$ and
[Fe/H]222these parameters come from a grid search, so that for small changes
in $\log{g}$ ($\Delta\log{g}<0.01$) the model with the best matching $T_{\rm
eff}$ and [Fe/H] will be the same. . The obtained parameters for each star are
presented in Table 2.

Table 2: Stellar propertiesa of HATS-38 and WASP-139.

Parameter | HATS-38 | WASP-139 | Reference
---|---|---|---
RA (J2015.5) | 10h17m05.05s | 03h18m14.92s | Gaia DR2
Dec (J2015.5) | -25h16m34.67s | -41d18m07.28s | Gaia DR2
pmRA (mas/yr) | -21.75$\pm$0.07 | -15.99$\pm$0.03 | Gaia DR2
pmDEC (mas/yr) | -7.54$\pm$0.07 | 24.51$\pm$0.05 | Gaia DR2
$\pi$ (mas) | 2.88$\pm$0.04 | 4.68$\pm$0.02 | Gaia DR2
T (mag) | 11.813$\pm$0.007 | 11.728$\pm$0.006 | TICv8
B (mag) | 12.76$\pm$0.42 | 13.330$\pm$0.016 | APASSb
V (mag) | 11.967$\pm$0.030 | 12.456$\pm$0.046 | APASS
G (mag) | 12.278$\pm$0.0002 | 12.271$\pm$0.0001 | Gaia DR2c
GBP (mag) | 12.649 $\pm$ 0.0002 | 12.733 $\pm$ 0.0002 | Gaia DR2
GRP (mag) | 11.761 $\pm$ 0.0002 | 11.677 $\pm$ 0.0002 | Gaia DR2
J (mag) | 11.184$\pm$0.026 | 10.982$\pm$0.023 | 2MASSd
H (mag) | 10.850$\pm$0.024 | 10.575$\pm$0.023 | 2MASS
Ks (mag) | 10.768$\pm$0.024 | 10.472$\pm$0.021 | 2MASS
$T_{\rm eff}$ (K) | 5662 $\pm$ 80 | 5233$\pm$60 | This work
$\log{g}$ (dex) | 4.34$\pm$0.02 | 4.56 $\pm$ 0.02 | This work
$[$Fe/H$]$ (dex) | 0.0 $\pm$ 0.05 | 0.03$\pm$0.05 | This work
$v\sin{i}$ (km/s) | 3.26$\pm$0.3 | 2.49 $\pm$0.3 | This work
$M_{\star}$ ($M_{\odot}$) | 0.92$\pm$0.03 | 0.84$\pm$0.03 | This work
$R_{\star}$ ($R_{\odot}$) | 1.08$\pm$0.02 | 0.810$\pm$0.008 | This work
L⋆ (L⊙) | 1.11$\pm$0.06 | 0.44$\pm$ 0.02 | This work
AV (mag) | 0.12$\pm$0.08 | 0.09$\pm$0.06 | This work
Age (Gyr) | 10.3${}_{-2.1}^{+1.6}$ | 6.2${}^{+3.1}_{-3.4}$ | This work
$\rho_{\star}$ (g/cm3) | 1.02${}_{-0.05}^{+0.08}$ | 2.25$\pm$0.13 | This work

Note. — a The stellar parameters computed in this work do not consider
possible systematic differences among different stellar evolutionary models
(Tayar et al., 2022) and have underestimated uncertainties, b Munari et al.
(2014),c Gaia Collaboration et al. (2018),d Skrutskie et al. (2006).

##  4 Photometric Analysis

In order to update the orbital ephemeris of HATS-38 b and WASP-139 b, and look
for Transit Timing Variations (TTVs), we performed a photometric analysis with
juliet (Espinoza et al., 2019) for each planet. juliet uses batman (Kreidberg,
2015) for the transit model and the dynesty dynamic nested sampler (Speagle,
2020) to perform bayesian analysis and explore the likelihood space to obtain
posterior probability distributions. We placed uninformative priors on the
transit parameters $R_{p}/R_{\star}$ and $b$, with an informative prior on the
stellar density, that was constrained in Section 3. We sampled the limb
darkening parameters using the quadratic $q_{1}$ and $q_{2}$ limb darkening
formalism from Kipping (2013a) with uniform priors. We placed Gaussian priors
for each transit mid-point based on the expected values calculated from the
orbital period and $t_{0}$ from Jordán et al. (2020) and Hellier et al. (2017)
for HATS-38 b and WASP-139 b, respectively, placing a large width of 0.1 days
on the Gaussian prior standard deviation to not impact the derived transit
midpoints. To account for variability and systematic noise in the TESS light
curves, we included a Matern-3/2 Gaussian Process (GP) as implemented in
celerite (Foreman-Mackey et al., 2017) and available in juliet. Each TESS year
(i.e., all combined sectors of each year) had its own GP kernel to account for
differences in variability captured in different epochs and cadences. For the
MOANA/ES light curve, we also included a linear trend using the airmass as
input as the observed flux was seen to be correlated with the airmass of the
observations. From this analysis, we ruled out the presence of TTVs greater
than $\sim 10$ minutes for HATS-38 b and $\sim 3$ minutes for WASP-139 b.
Additionally, we obtained detrended TESS and MOANA/ES light curves, and
ephemeris consistent with previous efforts that used different input data from
HATSouth (Bakos et al., 2013) and WASP (Pollacco et al., 2006).

##  5 Joint fit

To precisely constrain the parameters of HATS-38 b and WASP-139 b, and their
orbits, we broadly followed the methodology of Espinoza-Retamal et al. (2023),
which we have implemented in a publicly available code named
ironman333https://github.com/jiespinozar/ironman. Since this is the first time
we use the code, we describe it here in detail. ironman is a Python-based code
that allows the user to fit transit light curves, Keplerian RVs, and the RM
effect, all together in order to better constrain the parameters of a given
system. To model the transit light curves, ironman uses batman (Kreidberg,
2015). The default parameters required to model a light curve are the orbital
period ($P$), time of mid-transit ($t_{0}$), radius ratio ($R_{p}/R_{\star}$),
orbital inclination ($i$), scaled semimajor axis ($a/R_{\star}$), eccentricity
($e$), argument of periastron ($\omega$), limb darkening coefficients of the
instrument/band ($q_{1}^{\rm inst}$ and $q_{2}^{\rm inst}$), and the jitter
term ($\sigma_{\rm inst}$). In this case $q_{1}^{\rm inst}$ and $q_{2}^{\rm
inst}$ are the limb darkening parameters from Kipping (2013a). Additionally,
ironman can accept the impact parameter ($b$) and the stellar density
($\rho_{\star}$) instead of $i$ and $a/R_{\star}$, respectively. To model the
RVs, ironman uses rmfit (Stefansson et al., 2020, 2022), which uses radvel
(Fulton et al., 2018) to get the Keplerian orbit, and the framework from
Hirano et al. (2010) to model the RM effect. The required parameters to model
the Keplerian orbit of the planet are $P$, $t_{0}$, $e$, $\omega$, RV
semiamplitude ($K$), RV offset ($\gamma_{\rm inst}$), RV linear trend
($\dot{\gamma}$), RV quadratic trend ($\ddot{\gamma}$), and $\sigma_{\rm
inst}$. Finally, to model the RM effect, the sky-projected obliquity
($\lambda$), the projected rotational velocity of the star
($v\sin{i_{\star}}$), and the intrinsic linewidth $\beta_{\rm inst}$—which
accounts for instrumental and macroturbulence broadening—are also required.
Further, ironman can also accept a different parametrization of the RM effect,
sampling the RM model using the stellar inclination ($\cos{i_{\star}}$),
stellar rotational period ($P_{\rm rot}$), and stellar radius ($R_{\star}$)
instead of $v\sin{i_{\star}}$, in order to constrain the stellar inclination
and the true obliquity, following the parametrization used in Stefansson et
al. (2022). If this last parametrization is used, ironman estimates the true
3D obliquity $\psi$ as:

$\cos{\psi}=\cos{i_{\star}}\cos{i}+\sin{i_{\star}}\sin{i}\cos{\lambda}.$ (1)

To estimate the Bayesian posteriors, ironman uses the dynesty dynamic nested
sampler (Speagle, 2020). The number of live points to consider and threads to
be used can be specified by the user. ironman accepts uniform, log-uniform,
normal, and truncated normal prior distributions for the different parameters,
which can also be held fixed.

Figure 3: Sky-projected obliquities of Neptunes ($10<M_{p}/M_{\oplus}<50$ or
$2<R_{p}/R_{\oplus}<6$) as a function of stellar effective temperature (blue).
Neptunes on compact multiplanet systems are shown as squares, while isolated
Neptunes are shown as circles. Jupiters (not Neptunes with
$R_{p}/R_{\oplus}>6$) and super-Earths (not Neptunes with
$R_{p}/R_{\oplus}<2$) from TEPCat (Southworth, 2011) are shown in orange and
grey, respectively.

With ironman, we analyzed both systems, including all the observations
presented in Section 2. For this analysis, we only considered the detrended
TESS light curves (see Section 4) close to the transits, in windows of $\sim
6$ hours around the transit midpoint, to reduce the computational cost. To
account for the instrumentation offsets and systematics, we included
independent RV offsets and log-uniform jitter terms for each instrument. We
placed uninformative priors for almost all parameters, except for the orbital
periods and times of mid-transit, which had been constrained in Section 4, and
for the stellar densities that were constrained in Section 3. Likewise for
$\beta$ where we considered an instrumental broadening of 2.15 km/s because of
the ESPRESSO resolution, and we used the macroturbulence law from Valenti &
Fischer (2005), which yields a broadening of $\sim 4.15$ and $\sim 4.8$ km/s
for $T_{\rm eff}=5662$ and $5233$ K, respectively. We added the instrumental
and macroturbulence broadening in quadrature to set our priors for each target
respectively, with an uncertainty of 2 km/s. The resulting posteriors of our
fits for HATS-38 b and WASP-139 b are shown in Table 3.

Figures 1 and 2 show the different datasets together with the best-fit model.
We found that both planets lie in a nearly-polar orbit, with sky-projected
obliquities $\lambda=-108_{-16}^{+11}$ and $-85.6_{-4.2}^{+7.7}$ deg for
HATS-38 b and WASP-139 b, respectively. Also, we found that both planets are
consistent with having eccentric orbits. This was not noted in Jordán et al.
(2020) and Hellier et al. (2017), where they imposed circular orbits. We
tested fitting circular orbits as well, but the eccentric fit had a higher
evidence value, with a $\Delta\log{Z}>2$ in both cases, so it is the one we
adopted. We also tried fitting a long-term linear RV slope to the data, but it
resulted in being consistent with zero so we fixed it to that value. Overall,
the rest of the parameters are consistent with the ones reported in Jordán et
al. (2020) and Hellier et al. (2017).

Finally, we evaluated the possibility of a correlation between the measured
ESPRESSO RVs and stellar activity. Using actin2 (Gomes da Silva et al., 2018,
2021) and the ESPRESSO DAS pipeline (v. 1.3.7, Cupani et al., 2015) we
calculated a number of activity indexes including H$\alpha$, $R_{HK}$, $S$,
the bisector span, and the FWHM of the CCFs. We do not observe any significant
correlation between the RVs and any of the indexes in both datasets.

##  6 Discussion

In this work, we have characterized the orbits HATS-38 b and WASP-139 b, and
found that they have nearly polar sky-projected obliquities and residual non-
zero eccentricities, thus joining the growing sample of isolated (without
nearby companions within 4 times the orbital period) low-density Neptunes in
polar and eccentric orbits that include HAT-P-11 b (Sanchis-Ojeda & Winn,
2011), WASP-107 b (Dai & Winn, 2017; Rubenzahl et al., 2021), GJ 436 b
(Bourrier et al., 2018, 2022), and GJ 3470 b (Stefansson et al., 2022).

Table 3: Summary of posteriors of the joint fits.  Parameter | Description | HATS-38 b | WASP-139 b
---|---|---|---
$\lambda$ | Sky-projected obliquity (deg) | $-108_{-16}^{+11}$ | $-85.6_{-4.2}^{+7.7}$
$v\sin{i}$ | Projected rotational velocity (km/s) | $1.9_{-0.6}^{+1.1}$ | $1.8_{-0.8}^{+1.6}$
$P$ | Orbital period (days) | $4.37504_{-0.00002}^{+0.00002}$ | $5.9242705_{-0.0000007}^{+0.0000008}$
$t_{0}$ | Transit midpoint (BJD) | $2457786.4105_{-0.0005}^{+0.0005}$ | $2456870.9582_{-0.0002}^{+0.0002}$
$\rho_{\star}$ | Stellar density (g/cm3) | $1.02_{-0.08}^{+0.08}$ | $2.28_{-0.13}^{+0.13}$
$b$ | Impact parameter | $0.42_{-0.21}^{+0.14}$ | $0.19_{-0.09}^{+0.11}$
$R_{p}/R_{\star}$ | Radius ratio | $0.057_{-0.001}^{+0.001}$ | $0.099_{-0.001}^{+0.001}$
$e$ | Eccentricity | $0.112_{-0.070}^{+0.072}$ | $0.103_{-0.041}^{+0.050}$
$\omega$ | Argument of periastron (deg) | $159_{-56}^{+70}$ | $140_{-40}^{+20}$
$K$ | RV semiamplitude (m/s) | $8.6_{-1.8}^{+2.0}$ | $14.9_{-2.2}^{+2.1}$
$a/R_{\star}$ | Scaled semimajor axis | $10.1_{-0.3}^{+0.2}$ | $16.2_{-0.3}^{+0.3}$
$i$ | Orbital inclination (deg) | $87.52_{-0.59}^{+1.14}$ | $89.28_{-0.37}^{+0.35}$
$a$ | Semimajor axis (au) | $0.051_{-0.002}^{+0.002}$ | $0.061_{-0.001}^{+0.001}$
$R_{p}$ | Planet radius ($R_{\oplus}$) | $6.7_{-0.2}^{+0.2}$ | $8.8_{-0.1}^{+0.1}$
$M_{p}$ | Planet mass ($M_{\oplus}$) | $20.7_{-4.3}^{+4.8}$ | $37.6_{-5.5}^{+5.4}$
$\rho_{p}$ | Planet mean density (g/cm3) | $0.38_{-0.09}^{+0.10}$ | $0.31_{-0.05}^{+0.05}$
$q_{1}^{\rm ESPRESSO}$ | ESPRESSO linear limb darkening parameter | $0.74_{-0.28}^{+0.18}$ | $0.54_{-0.32}^{+0.29}$
$q_{2}^{\rm ESPRESSO}$ | ESPRESSO quadratic limb darkening parameter | $0.65_{-0.37}^{+0.25}$ | $0.43_{-0.31}^{+0.36}$
$\beta_{\rm ESPRESSO}$ | Intrinsic ESPRESSO stellar line width (km/s) | $4.4_{-2.0}^{+1.9}$ | $5.7_{-1.8}^{+1.8}$
$\gamma_{\rm ESPRESSO}$ | ESPRESSO RV offset (m/s) | $4139_{-1}^{+1}$ | $-13004_{-1}^{+1}$
$\gamma_{\rm HARPS}$ | HARPS RV offset (m/s) | $4144_{-1}^{+1}$ | $-13011_{-2}^{+1}$
$\gamma_{\rm PFS}$ | PFS RV offset (m/s) | $-1_{-2}^{+4}$ | -
$\gamma_{\rm FEROS}$ | FEROS RV offset (m/s) | $1_{-6}^{+5}$ | -
$\gamma_{\rm CORALIE1}$ | CORALIE1 RV offset (m/s) | - | $-13013_{-4}^{+4}$
$\gamma_{\rm CORALIE2}$ | CORALIE2 RV offset (m/s) | - | $-12988_{-8}^{+8}$
$\sigma_{\rm ESPRESSO}$ | ESPRESSO RV jitter (m/s) | $0.1_{-0.1}^{+0.5}$ | $0.1_{-0.1}^{+0.6}$
$\sigma_{\rm HARPS}$ | HARPS RV jitter (m/s) | $0.6_{-0.6}^{+3.9}$ | $0.2_{-0.2}^{+1.3}$
$\sigma_{\rm PFS}$ | PFS RV jitter (m/s) | $4.6_{-4.4}^{+5.1}$ | -
$\sigma_{\rm FEROS}$ | FEROS RV jitter (m/s) | $15.9_{-4.9}^{+6.1}$ | -
$\sigma_{\rm CORALIE1}$ | CORALIE1 RV jitter (m/s) | - | $0.4_{-0.4}^{+6.5}$
$\sigma_{\rm CORALIE2}$ | CORALIE2 RV jitter (m/s) | - | $0.2_{-0.2}^{+3.5}$
$q_{1}^{\rm TESS}$ | TESS linear limb darkening parameter | $0.36_{-0.18}^{+0.25}$ | $0.34_{-0.11}^{+0.14}$
$q_{2}^{\rm TESS}$ | TESS quadratic limb darkening parameter | $0.16_{-0.12}^{+0.25}$ | $0.29_{-0.12}^{+0.17}$
$q_{1}^{r^{\prime}}$ |  $r^{\prime}$ linear limb darkening parameter | $0.19_{-0.11}^{+0.19}$ | -
$q_{2}^{r^{\prime}}$ |  $r^{\prime}$ quadratic limb darkening parameter | $0.28_{-0.21}^{+0.37}$ | -
$q_{1}^{i^{\prime}}$ |  $i^{\prime}$ linear limb darkening parameter | $0.16_{-0.07}^{+0.11}$ | -
$q_{2}^{i^{\prime}}$ |  $i^{\prime}$ quadratic limb darkening parameter | $0.6_{-0.34}^{+0.27}$ | -
$q_{1}^{I+Z}$ |  $I+z$ linear limb darkening parameter | - | $0.73_{-0.18}^{+0.16}$
$q_{2}^{I+Z}$ |  $I+z$ quadratic limb darkening parameter | - | $0.08_{-0.05}^{+0.09}$
$q_{1}^{NGTS}$ |  $NGTS$ linear limb darkening parameter | - | $0.52_{-0.13}^{+0.14}$
$q_{2}^{NGTS}$ |  $NGTS$ quadratic limb darkening parameter | - | $0.44_{-0.11}^{+0.15}$
$\sigma^{\rm TESS}_{\rm Y1}$ | TESS Year 1 photometric jitter (ppm) | $19_{-16}^{+89}$ | $14_{-11}^{+74}$
$\sigma^{\rm TESS}_{\rm Y3}$ | TESS Year 3 photometric jitter (ppm) | $12_{-9}^{+66}$ | $11_{-9}^{+73}$
$\sigma^{\rm TESS}_{\rm Y5}$ | TESS Year 5 photometric jitter (ppm) | $16_{-13}^{+79}$ | -
$\sigma^{r^{\prime}}_{\rm MOANA/ES}$ | MOANA $r^{\prime}$ photometric jitter (ppm) | $3005_{-135}^{+147}$ | -
$\sigma^{r^{\prime}}_{\rm CHAT}$ | CHAT $r^{\prime}$ photometric jitter (ppm) | $898_{-59}^{+65}$ | -
$\sigma^{i^{\prime}}_{\rm LCO/SAAO}$ | LCO/SAAO $i^{\prime}$ photometric jitter (ppm) | $51_{-46}^{+215}$ | -
$\sigma^{i^{\prime}}_{\rm LCO/CTIO}$ | LCO/CTIO $i^{\prime}$ photometric jitter (ppm) | $16_{-13}^{+110}$ | -
$\sigma^{I+z}_{\rm TRAPPIST}$ | TRAPPIST 2014 $I+z$ photometric jitter (ppm) | - | $111_{-105}^{+705}$
$\sigma^{I+z}_{\rm TRAPPIST}$ | TRAPPIST 2015 $I+z$ photometric jitter (ppm) | - | $2320_{-138}^{+142}$
$\sigma^{NGTS}_{\rm EulerCAM}$ | EulerCAM $NGTS$ photometric jitter (ppm) | - | $593_{-87}^{+86}$

In Figure 3, we compare our results for HATS-38 b and WASP-139 b to the sample
of Neptunes with measured sky-projected obliquity from
TEPCat444https://www.astro.keele.ac.uk/jkt/tepcat/ (Southworth, 2011),
plotting the sky-projected obliquity versus stellar effective temperature. In
this plot, we can see that Neptunes in compact multiplanet systems are all
consistent with being well-aligned, except for the notable case of HD 3167c
(Dalal et al., 2019; Bourrier et al., 2021). In turn, isolated Neptunes can be
either polar or well-aligned. This dichotomy has recently been discussed by
Radzom et al. (2024).

Next, we discuss some of the properties of HATS-38 b and WASP-139 b as well as
the overall distribution of stellar obliquities of the close-in Neptunes, some
of the possible origins of these planets, and put some constraints on the
possibility of having companions.

###  6.1 A possible preponderance of polar Neptunes?

To better quantify the possible preponderance of polar Neptunes in the sample,
using the two new measurements we have provided here, we used a hierarchical
Bayesian model (e.g., Hogg et al., 2010) for the underlying obliquity
distribution (e.g., Morton & Winn, 2014; Muñoz & Perets, 2018). In particular,
we used the framework555https://github.com/jiayindong/obliquity presented by
Dong & Foreman-Mackey (2023), which models the underlying distribution of
$\cos{\psi}$ across an exoplanet population using a mixture model of two Beta
distributions (e.g., Gelman et al., 2014). This hierarchical Bayesian modeling
framework allows the user to derive the true 3D obliquity distribution from
observed sky-projected obliquities and has the capacity to capture anything
from an isotropic distribution to a strongly bimodal population.

Both Beta distribution components are modeled using 3 parameters each: $w$,
$\mu$, and $\kappa$. The parameter $w$ describes the weight of the component,
while $\mu$ and $1/\kappa$ correspond to the mean and variance of each beta
distribution component. The greater the value of $\kappa$, the smaller the
variance, i.e., the distribution is more concentrated. The $\mu$ and $\kappa$
parameters can be related to the typical $\alpha$ and $\beta$ parameters of a
beta distribution: $\alpha=\mu\kappa$ and $\beta=(1-\mu)\kappa$.

Figure 4: Inferred stellar obliquity distribution of Neptunes (blue) and
Jupiters (orange) with sky-projected stellar obliquity measurements. This
inference was done following the methodology from Dong & Foreman-Mackey (2023)
using the observed sky-projected obliquities available in the
literature—combined with the two sky-projected obliquity measurements
presented here. The lines and shaded regions represent the median and
$1\sigma$ uncertainties of the inferred distributions. We show the results for
a) a $\mathcal{N}(0,3)$ prior on the $\log{\kappa}$ parameter and b) a
$\mathcal{U}(-4,10)$ prior on $\log{\kappa}$.

To derive the underlying $\cos{\psi}$ distribution, we ran two different
hierarchical models with different priors on the compactness parameter
$\kappa$: first, a model using informative priors on $\kappa$ following Dong &
Foreman-Mackey (2023), and second, a fit using uninformative priors on
$\kappa$. For the first model, we followed Dong & Foreman-Mackey (2023), and
assumed an informative normal prior on $\log{\kappa}$ with a mean of 0 and a
standard deviation of 3. Second, as an additional test—given the small size of
the sample—we ran another model with an uninformative uniform prior between -4
and 10 for $\log{\kappa}$. In both cases, we used uniform priors on the
location parameter $\mu$. The results of both of these model fits are shown in
Table 4, and in Figure 4.

We have found tentative evidence for a dichotomy in the orbital alignment of
the sample of 27 Neptunes. Their orbits are preferentially well-aligned or
polar. The relative weight of the two peaks depends on the priors we assumed,
being $w_{0}=0.44_{-0.10}^{+0.11}$ and $w_{1}=0.56_{-0.11}^{+0.10}$ in the
informative case, and $w_{0}=0.39_{-0.31}^{+0.24}$ and
$w_{1}=0.61_{-0.24}^{+0.31}$ in the uninformative case. These weights suggest
that about half of the Neptunes are expected to be in polar orbits while the
other half is expected to be in well-aligned orbits. This is also reflected in
the as-measured sky-projected obliquity sample: with the addition of our two
new polar measurements, 13 Neptunes have $|\lambda|>50$ deg, and the remaining
14 have $|\lambda|<50$ deg. For both sets of priors, the polar peak is around
$\cos{\psi}\sim-0.18$ while the aligned peak is around $\cos{\psi}\sim 1$.
However, we see that the concentration of the polar peak and its uncertainty
are strongly dependent on the priors used on $\log\kappa$, where the two
different runs yield $\log{\kappa_{0}}=3.08_{-1.47}^{+1.64}$ in the
informative case and $\log{\kappa_{0}}=4.61_{-3.63}^{+5.15}$ in the
uninformative case. For reference, $\log{\kappa_{0}}$ values of $\lesssim 1.0$
correspond to isotropic distributions as seen from the $\log{\kappa_{0}}$
values of $1.24_{-0.92}^{+1.11}$ and $0.58_{-0.91}^{+3.26}$ for the Jupiter
population in the informative and non-informative cases, respectively (see
Table 4). If we define a $\log(\kappa_{0})=1$ to represent an isotropic
distribution, we see that the Neptunes are $\sim$1.5$\sigma$ non-isotropic in
the informative case, and consistent with isotropic at $\sim$1$\sigma$ in the
non-informative case.

In contrast, the $\cos{\psi}$ distribution of Jupiters (171 systems,
representing the entire sample with no cut in $T_{\rm eff}$) is not sensitive
to the prior assumptions. Figure 4 shows that the distribution has a peak at 0
deg representing the well-aligned systems, with an almost isotropic tail and
no significative clustering at 90 deg. However, the respective weights of each
beta component are different for the different priors. In the informative
case, we have $w_{0}=0.33_{-0.06}^{+0.07}$ and $w_{1}=0.67_{-0.07}^{+0.06}$,
suggesting that $\sim 1/3$ of the hot Jupiters are in misaligned orbits and
the other $\sim 2/3$ are in well-aligned orbits. These results are consistent
with the values derived by Dong & Foreman-Mackey (2023) for the sample of
exoplanets from Albrecht et al. (2022), which is dominated by hot Jupiters.
However, with the uninformative priors, the weights are less constrained with
$w_{0}=0.41\pm 0.40$ and $w_{1}=0.59\pm 0.40$ maybe suggesting that we could
model the distribution using only 1 component.

The observed difference in the obliquity distribution for the sample of
Neptunes and Jupiters—with the former showing a possible preference for a
polar + aligned obliquity dichotomy, while the latter shows a preference for
an aligned + isotropic dichotomy—likely suggests that hot Jupiters and hot
Neptunes could have formed through different formation channels. However, the
sample of Neptunes is small compared to the one of Jupiters (27 versus 171
systems), and the exact parameters of the bimodal distribution are different
with different priors. According to Dong & Foreman-Mackey (2023), in their
tests, prior distribution choices did not strongly affect the derived
distributions for samples of more than $\sim 50$ systems. Therefore, more
measurements of the stellar obliquity of Neptunes will be required to further
confirm or rule out these results. Future efforts like this one, and programs
like The Desert-Rim Exoplanets Atmosphere and Migration (DREAM, Bourrier et
al., 2023), will be crucial for increasing the size of the sample and
understanding the formation and evolution of the hot Neptunes.

###  6.2 Companions

The baseline of the RV measurements used in this work for HATS-38 and WASP-139
are $\sim 7$ and $\sim 15$ years respectively, and both residuals show no
clear evidence of long-period companions (see panel c in Figures 1 and 2). In
order to estimate what kind of planetary companions can be ruled out by the
current RVs, we performed population synthesis injection tests. For a
synthetic population of possible companions, we estimated their RV signals and
checked if they were consistent with the residuals of the RV model for HATS-38
b and WASP-139 b.

In both cases, we considered planetary companions with masses between 0.3 and
50 $M_{J}$. The eccentricity of the planets was drawn from a uniform
distribution between 0.2 and 0.5, considering that cold Jupiters tend to be
moderately eccentric (see, e.g., Kipping, 2013b). The cosine of the orbital
inclinations of the companions was drawn uniformly between 0 and 1666Similar
results are obtained with the extreme case where all companions are edge-on..
Finally, the semimajor axes were drawn uniformly from 0.5 and 100 au. With
those distributions, we estimated the RV signals of different companions.
Figure 5 shows regions in the mass-semimajor axis plane ruled out at different
confidence levels by the RV residuals. In both systems, the RVs rule out the
existence of any planetary companions with masses between 0.3 and 50 $M_{J}$
within $\sim 10$ au at $5\sigma$ confidence level.

Additionally to these possible planetary companions, we note that both stars
appear to have stellar companions. Jordán et al. (2020) reported a possible
M-dwarf companion to HATS-38, with a mass of 0.1 $M_{\odot}$ and a current
projected separation of $\sim 2100$ au. However, the source is not included in
the Gaia catalog, given its low expected G magnitude of 23. So it can also be
an extragalactic source, an earlier M-dwarf star that is in the background of
HATS-38, or a foreground brown dwarf. As for WASP-139, El-Badry et al. (2021)
reported a comoving M-dwarf companion with a mass of 0.1 $M_{\odot}$ at a
current projected separation of $\sim 6400$ au. In both cases, the M-dwarf
companion could have influenced the architecture of the system, as discussed
in the next Section.

Table 4: Summary of posteriors for the hyperparameters of the derived two-
component $\cos{\psi}$ distribution.

|  | Misaligned Component | Well-aligned Component |
---|---|---|---|---
Sample | Systems | $w_{0}$ | $\mu_{0}$ | $\log{\kappa_{0}}$ | $w_{1}$ | $\mu_{1}$ | $\log{\kappa_{1}}$ | Reference
Prior on $\log\kappa$ of $\mathcal{N}(0,3)$a |  |  |  |  |  |  |  |
Neptunes | 27 | $0.44_{-0.10}^{+0.11}$ | $0.41\pm 0.05$ | $3.08_{-1.47}^{+1.64}$ | $0.56_{-0.11}^{+0.10}$ | $0.97_{-0.03}^{+0.01}$ | $3.01_{-2.17}^{+1.78}$ | This work
Jupiters | 171 | $0.33_{-0.06}^{+0.07}$ | $0.47_{-0.07}^{+0.08}$ | $1.24_{-0.92}^{+1.11}$ | $0.67_{-0.07}^{+0.06}$ | $0.98_{-0.01}^{+0.01}$ | $3.01_{-1.80}^{+1.67}$ | This work
A2022b | 161 | $0.281\pm 0.085$ | $0.434\pm 0.088$ | $1.44\pm 1.33$ | $0.719\pm 0.085$ | $0.976\pm 0.022$ | $2.65\pm 1.88$ | D&FM2023c
Prior on $\log\kappa$ of $\mathcal{U}(-4,10)$a |  |  |  |  |  |  |  |
Neptunes | 27 | $0.39_{-0.31}^{+0.24}$ | $0.40_{-0.14}^{+0.12}$ | $4.61_{-3.63}^{+5.15}$ | $0.61_{-0.24}^{+0.31}$ | $0.94_{-0.23}^{+0.05}$ | $1.61_{-2.08}^{+5.61}$ | This work
Jupiters | 171 | $0.41\pm 0.40$ | $0.56_{-0.38}^{+0.22}$ | $0.58_{-0.91}^{+3.26}$ | $0.59\pm 0.40$ | $0.99_{-0.19}^{+0.01}$ | $3.88_{-3.99}^{+5.13}$ | This work

Note. — a $\mathcal{N}(\mu,\sigma)$ denotes a normal prior with median $\mu$
and standard deviation $\sigma$ while $\mathcal{U}(a,b)$ denotes a uniform
prior between $a$ and $b$, bAlbrecht et al. (2022), cDong & Foreman-Mackey
(2023).

###  6.3 Possible Origin of the Polar Orbits

We discuss some of the proposed explanations for the origin of the polar
orbits of hot Neptunes and whether these are able to explain not only the
stellar obliquities of HATS-38 b and WASP-139 b, but also their residual
eccentricities and constraints on long-period companions. We argue that
primordial misalignments and in situ excitation are ruled out by these
observational constraints, while high-eccentricity migration from $\gtrsim 2$
au remains as a possible formation pathway.

#### Primordial disk misalignment

Highly misaligned orbits can be produced through a primordial tilt of the
protoplanetary disk by gravitational torques from a binary stellar companion
(e.g., Batygin, 2012; Lai, 2014; Zanazzi & Lai, 2018). This scenario seems
possible given that both stars appear to have distant M-dwarf companions at
projected separations of $\sim 2100$ au and $\sim 6400$ au (Section 6.2).
However, at these large separations a precession cycle of the disk is
long777For reference, a protoplanetary disk extending out to 200 au would
complete a precession cycle in $\sim 15$ Myr due to a M-dwarf companion with
semi-major axis of 2,000 au (e.g., Equation (12) in Zanazzi & Lai 2018).
compared to the disk dispersal timescale to resonantly excite the stellar
obliquities (Zanazzi & Lai, 2018). Alternatively, the stellar obliquity may be
excited by the disk through magnetic disk torquing in a young accreting star
(Lai et al., 2011). Although these may still be plausible processes to
misalign the parent disks of HATS-38 b and WASP-139 b to the equator of their
host stars, they are unable to explain the observed eccentricities. Therefore,
we disfavor primordial misalignments as an explanation.

Figure 5: Mass versus semimajor axis diagram for the possible outer companions
to a) HATS-38 b and b) WASP-139 b. Different colors indicate regions where we
can rule out the companion at different $\sigma$ levels based on the RV
residuals. Yellow dotted-dashed lines show $\eta_{\rm GR}=3$ (Equation 4),
demarking the regions where the companion can excite the eccentricities
secularly. Cyan dotted-dashed lines show $\eta_{\star}=1$ (Equation 5) for
$J_{2}=10^{-7}$ (corresponding to rotational periods of $\sim 20-30$ days),
demarking the regions where the companion can secularly excite obliquities.
Blue dashed lines show $\tau_{\rm ZLK}=1$ Myr (Equation 3), demarking the
regions where the companion assisted by a dispersing protoplanetary disk can
resonantly excite obliquities (adiabatic resonance capture). Black dashed
lines show $\eta_{\rm GR}=1/2$, demarking the regions where the companion can
induce high-eccentricity migration if it started at $a_{\rm b,0}=2$ au. For
any of the mechanisms discussed above to play a role, the companion needs to
reside above the respective line; for companions below the line, the
respective mechanism is excluded. From both figures, we see that the current
RV observations are only compatible with the high-eccentricity migration
scenario.

#### In-situ planet-planet scattering

The excited orbital state of planet b could have resulted from gravitational
instabilities near its current orbital separation. As no other planets (with
similar or large masses) are detected, such planets should have been ejected,
leaving planet b in an eccentric/inclined orbit. We quantify the likelihood of
this scenario by the Safronov-like parameter $\theta$ defined as the ratio of
the escape velocity of the planet ($V_{\rm esc}$) and the orbital velocity
($V_{\rm orb}$):

$\theta^{2}\equiv\frac{V_{\rm esc}^{2}}{V_{\rm
orb}^{2}}=\frac{2M_{p}}{M_{\star}}\cdot\frac{a}{R_{p}}.$ (2)

For HATS-38 b and WASP-139 b, we obtain $\theta^{2}\simeq 0.02$ and
$\theta^{2}\simeq 0.04$, respectively, implying in both cases that scattering
leads mainly to planetary collisions, not ejections. Thus, the excitation of
eccentricities and inclinations becomes unlikely at the planet’s orbital
separation (Petrovich et al., 2014).

#### In-situ excitation by a distant companion

A distant giant planet companion (planet "c") can secularly tilt the planetary
orbit and drive eccentricity growth on a time scale given by

$\tau_{\rm ZLK}=\frac{2P_{\rm b}^{2}}{3\pi P_{\rm
c}}\frac{m_{c}}{M_{\star}}\left(\frac{a_{c}}{a_{b}}\right)^{3},$ (3)

where ZLK stands for von Zeipel-Lidov-Kozai (von Zeipel, 1910; Lidov, 1962;
Kozai, 1962). The eccentricity growth can be quenched by relativistic apsidal
precession (rate defined by $\dot{\omega}_{\rm GR}$) and its strength relative
to the two-planet interaction can be quantified by the dimensionless parameter
as (e.g., Fabrycky & Tremaine, 2007; Liu et al., 2015)

$\eta_{\rm GR}=\dot{\omega}_{\rm GR}\tau_{\rm
ZLK}=\frac{4GM_{\star}}{c^{2}}\frac{a_{c}^{3}(1-e_{c}^{2})^{3/2}}{a_{b}^{4}}\frac{M_{\star}}{m_{c}}.$
(4)

For $\eta_{\rm GR}>3$ eccentricity excitation is quenched. In Figure 5, we
show that the allowed regions for HATS-38 b and WASP-139 b (pink dot-dashed
lines) are largely excluded by long-term RV measurements.

The inclinations can, in turn, be secularly excited by an inclined companion
and compete with the quadrupolar gravitational field with the host star
parametrized by $J_{2}$. We define the relative stellar quadrupole with
respect to the two-planet interaction as in Petrovich et al. (2020)

$\eta_{\star}\equiv\frac{2J_{2}M_{\star}}{m_{c}}\frac{R_{\star}^{2}a_{c}^{3}(1-e_{c}^{2})^{3/2}}{a_{b}^{5}}.$
(5)

For $\eta_{\star}>1$ a distant planet "c" is unable to drive the growth of the
stellar obliquity. In both panels of Figure 5 we plot $\eta_{\star}=1$
assuming $J_{2}=10^{-7}$ and the observed stellar masses and radii, and find
that we can rule out the presence of a planetary companion in most of the
regions where the necessary condition is satisfied.

A final possibility we assess is that the secular excitation occurred early
on, assisted by the dispersal of the birth protoplanetary disk, as proposed by
Petrovich et al. (2020). Here, the disk promotes the capture into an
inclination secular resonance (with possible eccentricity excitation) if the
dispersal of the disk is longer than the adiabatic timescale, which is defined
as $\tau_{\rm adia}=\tau_{\rm ZLK}(1+\eta_{\star})^{1/3}/I_{c,0}^{4/3}$ where
$I_{c,0}$ is the initial inclination of planet c relative to the disk. Since
$I_{c,0}\ll 1$, we can comfortably assume that $\tau_{\rm adia}\gg\tau_{\rm
ZLK}$ and use this bound to compare it with a typical dispersal timescale of
the disk of $1$ Myr. In Figure 5 we show the $\tau_{\rm ZLK}=1$ Myr, which
again, the condition $\tau_{\rm ZLK}<1$ Myr is largely ruled out by the radial
velocity constraints.

In summary, the secular excitation of the planet’s stellar obliquity at its
current orbital separation is largely disfavored by our RV constraints.

#### High-eccentricity migration driven by a distant companion

We have ruled out various mechanisms that can explain the excited orbital
states of the Neptunes at their current orbital distances. We turn to a final
possibility, namely excitation of their stellar obliquities and eccentricities
at wider orbital separations, followed by inward tidal migration.

High-eccentricity tidal migration driven by planet-planet scattering
accompanied by secular interactions (e.g., Nagasawa et al., 2008) has been
proposed as an explanation for the eccentric and misaligned orbits of hot
Neptunes (see, e.g., Bourrier et al., 2018; Correia et al., 2020; Stefansson
et al., 2022), including recent applications to the systems HAT-P-11 (Lu et
al., 2024) and WASP-107 (Yu & Dai, 2024). But, where should the migration have
started to still be consistent with the excluded regions by our RV
measurements?

A planet b can acquire an extreme eccentricity to allow for migration only if
$\eta_{\rm GR}\lesssim 1$. For reference, a maximum eccentricity of $e_{\rm
max}\gtrsim 0.95$ can be achieved only if $\eta_{\rm GR}>1/2$. In Figure 5 we
plot $\eta_{\rm GR}=1/2$ for an initial semimajor axis of the transiting
planet of $a_{\rm b,0}=2$ au, which coincides with regions excluded at only
$\gtrsim 1\sigma$ for companions with masses $\lesssim 3$ $M_{J}$. Given the
steep dependence with the initial semi-major axis ($\eta_{\rm GR}\propto
a_{\rm b,0}^{-4}$), we can conclude that ZLK migration could still be
consistent with excluded regions if $a_{\rm b,0}\gtrsim 2$ au.

If no companions are found, we are left with the wide-orbit binary companions
as the driver of ZLK migration (e.g., Fabrycky & Tremaine 2007; Naoz et al.
2012; Petrovich 2015), though migration should have started from an even wider
orbit to not be quenched by GR ($\eta_{\rm GR}<1/2$ for $a_{\rm b,0}\gtrsim 9$
au with $m_{c}=0.1$ $M_{\odot}$ and $a_{c}=2000$ au, see Section 6.2).
Additionally, the migrating planets should lack any other planetary companions
to avoid quenching the ZLK cycles, a state that could be achieved by planet-
planet scattering.

##  7 Conclusions and Summary

In this work, we have presented new observations of the hot Neptunes HATS-38 b
and WASP-139 b. We have jointly analyzed ESPRESSO observations of the RM
effect, produced by the transiting planets, with photometry and out-of-transit
radial velocities using the publicly available code ironman. We have concluded
that:

  * •

Both planets have nearly polar and eccentric orbits. We have measured a sky-
projected obliquity $\lambda=-108_{-16}^{+11}$ and $-85.6_{-4.2}^{+7.7}$ deg,
and an eccentricity $e=0.112_{-0.070}^{+0.072}$ and $0.103_{-0.041}^{+0.050}$
for HATS-38 b and WASP-139 b, respectively.

  * •

Neither systems show clear evidence of long-period companions in long-term RV
observations. Based on the residual RVs we have ruled out the presence of
planetary companion with masses $\sim 0.3-50$ $M_{J}$ within $\sim 10$ au at
$5\sigma$ confidence level.

  * •

Both planets join a growing population of isolated low-density Neptunes in
polar and eccentric orbits, which includes HAT-P-11 b (Sanchis-Ojeda & Winn,
2011), GJ 436 b (Bourrier et al., 2018, 2022), WASP-107 b (Dai & Winn, 2017;
Rubenzahl et al., 2021), and GJ 3470 b (Stefansson et al., 2022).

  * •

A possible explanation for the polar and eccentric architectures of these
systems is high-eccentricity migration produced by secular interactions from
yet undetected outer planetary companions or the observed wide-orbit stars. In
contrast, we disfavor in-situ excitation of inclinations and eccentricities,
as this would require companions at orbital distances we have ruled out with
out-of-transit RVs.

  * •

Through hierarchical Bayesian modeling of sky-projected obliquity
measurements, we found suggestive evidence for a preponderance of polar orbits
of hot Neptunes. In contrast, and similar to recent work by Siegel et al.
(2023) and Dong & Foreman-Mackey (2023), we do not observe significant
clustering of hot Jupiters in polar orbits. However, we note that the exact
obliquity distribution of Neptunes is sensitive to the choice of priors,
highlighting the need for additional obliquity measurements of Neptunes to
robustly compare the hot Neptune obliquity distribution to hot Jupiters.

## Acknowledgments

We would like to thank Coel Hellier for the TRAPPIST and EulerCAM light curves
of WASP-139. We also thank Jiayin Dong for her help with the implementation of
the Bayesian framework to derive the true 3D obliquity distributions. J.I.E.R
and C.P. acknowledge support from ANID BASAL project FB210003. J.I.E.R.
acknowledges support from ANID Doctorado Nacional grant 2021-21212378. C.P.
acknowledges support from FONDECYT Project 1210425, CASSACA grant CCJRF2105,
and ANID+REC Convocatoria Nacional subvención a la instalación en la Academia
convocatoria 2020 PAI77200076. A.J. and R.B. acknowledge support from ANID –
Millennium Science Initiative – ICN12_009. R.B. acknowledges support from
FONDECYT Project 1241963. A.J. acknowledges support from FONDECYT project
1210718. J.P.L. acknowledges co-funding by the European Union (ERC, FIERCE,
101052347), by FCT - Fundação para a Ciência e a Tecnologia - through national
funds, and by FEDER through COMPETE2020 - Programa Operacional Competitividade
e Internacionalização grants UIDB/04434/2020 and UIDP/04434/2020.

## References

  * Albrecht et al. (2012) Albrecht, S., Winn, J. N., Johnson, J. A., et al. 2012, ApJ, 757, 18, doi: 10.1088/0004-637X/757/1/18
  * Albrecht et al. (2022) Albrecht, S. H., Dawson, R. I., & Winn, J. N. 2022, PASP, 134, 082001, doi: 10.1088/1538-3873/ac6c09
  * Albrecht et al. (2021) Albrecht, S. H., Marcussen, M. L., Winn, J. N., Dawson, R. I., & Knudstrup, E. 2021, ApJ, 916, L1, doi: 10.3847/2041-8213/ac0f03
  * Bakos et al. (2013) Bakos, G. Á., Csubry, Z., Penev, K., et al. 2013, PASP, 125, 154, doi: 10.1086/669529
  * Batalha et al. (2013) Batalha, N. M., Rowe, J. F., Bryson, S. T., et al. 2013, ApJS, 204, 24, doi: 10.1088/0067-0049/204/2/24
  * Batygin (2012) Batygin, K. 2012, Nature, 491, 418, doi: 10.1038/nature11560
  * Bourrier et al. (2018) Bourrier, V., Lovis, C., Beust, H., et al. 2018, Nature, 553, 477, doi: 10.1038/nature24677
  * Bourrier et al. (2021) Bourrier, V., Lovis, C., Cretignier, M., et al. 2021, A&A, 654, A152, doi: 10.1051/0004-6361/202141527
  * Bourrier et al. (2022) Bourrier, V., Zapatero Osorio, M. R., Allart, R., et al. 2022, A&A, 663, A160, doi: 10.1051/0004-6361/202142559
  * Bourrier et al. (2023) Bourrier, V., Attia, O., Mallonn, M., et al. 2023, A&A, 669, A63, doi: 10.1051/0004-6361/202245004
  * Brahm et al. (2017a) Brahm, R., Jordán, A., & Espinoza, N. 2017a, PASP, 129, 034002, doi: 10.1088/1538-3873/aa5455
  * Brahm et al. (2017b) Brahm, R., Jordán, A., Hartman, J., & Bakos, G. 2017b, MNRAS, 467, 971, doi: 10.1093/mnras/stx144
  * Brahm et al. (2019) Brahm, R., Espinoza, N., Jordán, A., et al. 2019, AJ, 158, 45, doi: 10.3847/1538-3881/ab279a
  * Bressan et al. (2012) Bressan, A., Marigo, P., Girardi, L., et al. 2012, MNRAS, 427, 127, doi: 10.1111/j.1365-2966.2012.21948.x
  * Correia et al. (2020) Correia, A. C. M., Bourrier, V., & Delisle, J. B. 2020, A&A, 635, A37, doi: 10.1051/0004-6361/201936967
  * Cupani et al. (2015) Cupani, G., D’Odorico, V., Cristiani, S., et al. 2015, in Astronomical Society of the Pacific Conference Series, Vol. 495, Astronomical Data Analysis Software an Systems XXIV (ADASS XXIV), ed. A. R. Taylor & E. Rosolowsky, 289, doi: 10.48550/arXiv.1509.04901
  * Dai & Winn (2017) Dai, F., & Winn, J. N. 2017, AJ, 153, 205, doi: 10.3847/1538-3881/aa65d1
  * Dalal et al. (2019) Dalal, S., Hébrard, G., Lecavelier des Étangs, A., et al. 2019, A&A, 631, A28, doi: 10.1051/0004-6361/201935944
  * Dong & Foreman-Mackey (2023) Dong, J., & Foreman-Mackey, D. 2023, AJ, 166, 112, doi: 10.3847/1538-3881/ace105
  * El-Badry et al. (2021) El-Badry, K., Rix, H.-W., & Heintz, T. M. 2021, MNRAS, 506, 2269, doi: 10.1093/mnras/stab323
  * Espinoza et al. (2019) Espinoza, N., Kossakowski, D., & Brahm, R. 2019, MNRAS, 490, 2262, doi: 10.1093/mnras/stz2688
  * Espinoza-Retamal et al. (2023) Espinoza-Retamal, J. I., Brahm, R., Petrovich, C., et al. 2023, ApJ, 958, L20, doi: 10.3847/2041-8213/ad096d
  * Fabrycky & Tremaine (2007) Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298, doi: 10.1086/521702
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., & Angus, R. 2017, AJ, 154, 220, doi: 10.3847/1538-3881/aa9332
  * Freudling et al. (2013) Freudling, W., Romaniello, M., Bramich, D. M., et al. 2013, A&A, 559, A96, doi: 10.1051/0004-6361/201322494
  * Fulton et al. (2018) Fulton, B. J., Petigura, E. A., Blunt, S., & Sinukoff, E. 2018, PASP, 130, 044504, doi: 10.1088/1538-3873/aaaaa8
  * Gaia Collaboration et al. (2018) Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018, A&A, 616, A1, doi: 10.1051/0004-6361/201833051
  * Gelman et al. (2014) Gelman, A., Carlin, J. B., Stern, H. S., et al. 2014, Bayesian Data Analysis
  * Gomes da Silva et al. (2018) Gomes da Silva, J., Figueira, P., Santos, N., & Faria, J. 2018, The Journal of Open Source Software, 3, 667, doi: 10.21105/joss.00667
  * Gomes da Silva et al. (2021) Gomes da Silva, J., Santos, N. C., Adibekyan, V., et al. 2021, A&A, 646, A77, doi: 10.1051/0004-6361/202039765
  * Hellier et al. (2017) Hellier, C., Anderson, D. R., Collier Cameron, A., et al. 2017, MNRAS, 465, 3693, doi: 10.1093/mnras/stw3005
  * Hirano et al. (2014) Hirano, T., Sanchis-Ojeda, R., Takeda, Y., et al. 2014, ApJ, 783, 9, doi: 10.1088/0004-637X/783/1/9
  * Hirano et al. (2010) Hirano, T., Suto, Y., Taruya, A., et al. 2010, ApJ, 709, 458, doi: 10.1088/0004-637X/709/1/458
  * Hogg et al. (2010) Hogg, D. W., Myers, A. D., & Bovy, J. 2010, ApJ, 725, 2166, doi: 10.1088/0004-637X/725/2/2166
  * Jenkins et al. (2016) Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9913, Software and Cyberinfrastructure for Astronomy IV, ed. G. Chiozzi & J. C. Guzman, 99133E, doi: 10.1117/12.2233418
  * Jordán et al. (2020) Jordán, A., Bakos, G. Á., Bayliss, D., et al. 2020, AJ, 160, 222, doi: 10.3847/1538-3881/aba530
  * Kipping (2013a) Kipping, D. M. 2013a, MNRAS, 435, 2152, doi: 10.1093/mnras/stt1435
  * Kipping (2013b) —. 2013b, MNRAS, 434, L51, doi: 10.1093/mnrasl/slt075
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591, doi: 10.1086/108790
  * Kreidberg (2015) Kreidberg, L. 2015, PASP, 127, 1161, doi: 10.1086/683602
  * Lai (2014) Lai, D. 2014, MNRAS, 440, 3532, doi: 10.1093/mnras/stu485
  * Lai et al. (2011) Lai, D., Foucart, F., & Lin, D. N. C. 2011, MNRAS, 412, 2790, doi: 10.1111/j.1365-2966.2010.18127.x
  * Lidov (1962) Lidov, M. L. 1962, Planet. Space Sci., 9, 719, doi: 10.1016/0032-0633(62)90129-0
  * Lightkurve Collaboration et al. (2018) Lightkurve Collaboration, Cardoso, J. V. d. M., Hedges, C., et al. 2018, Lightkurve: Kepler and TESS time series analysis in Python, Astrophysics Source Code Library.  http://ascl.net/1812.013
  * Liu et al. (2015) Liu, B., Muñoz, D. J., & Lai, D. 2015, MNRAS, 447, 747, doi: 10.1093/mnras/stu2396
  * Lu et al. (2024) Lu, T., An, Q., Li, G., et al. 2024, arXiv e-prints, arXiv:2405.19511.  https://arxiv.org/abs/2405.19511
  * Masuda & Winn (2020) Masuda, K., & Winn, J. N. 2020, AJ, 159, 81, doi: 10.3847/1538-3881/ab65be
  * Mazeh et al. (2016) Mazeh, T., Holczer, T., & Faigler, S. 2016, A&A, 589, A75, doi: 10.1051/0004-6361/201528065
  * Modigliani et al. (2020) Modigliani, A., Freudling, W., Anderson, R. I., et al. 2020, in Astronomical Society of the Pacific Conference Series, Vol. 527, Astronomical Data Analysis Software and Systems XXIX, ed. R. Pizzo, E. R. Deul, J. D. Mol, J. de Plaa, & H. Verkouter, 667
  * Morton & Winn (2014) Morton, T. D., & Winn, J. N. 2014, ApJ, 796, 47, doi: 10.1088/0004-637X/796/1/47
  * Muñoz & Perets (2018) Muñoz, D. J., & Perets, H. B. 2018, AJ, 156, 253, doi: 10.3847/1538-3881/aae7d0
  * Munari et al. (2014) Munari, U., Henden, A., Frigo, A., et al. 2014, AJ, 148, 81, doi: 10.1088/0004-6256/148/5/81
  * Nagasawa et al. (2008) Nagasawa, M., Ida, S., & Bessho, T. 2008, ApJ, 678, 498, doi: 10.1086/529369
  * Naoz et al. (2012) Naoz, S., Farr, W. M., & Rasio, F. A. 2012, ApJ, 754, L36, doi: 10.1088/2041-8205/754/2/L36
  * Pepe et al. (2021) Pepe, F., Cristiani, S., Rebolo, R., et al. 2021, A&A, 645, A96, doi: 10.1051/0004-6361/202038306
  * Petrovich (2015) Petrovich, C. 2015, ApJ, 799, 27, doi: 10.1088/0004-637X/799/1/27
  * Petrovich et al. (2020) Petrovich, C., Muñoz, D. J., Kratter, K. M., & Malhotra, R. 2020, ApJ, 902, L5, doi: 10.3847/2041-8213/abb952
  * Petrovich et al. (2014) Petrovich, C., Tremaine, S., & Rafikov, R. 2014, ApJ, 786, 101, doi: 10.1088/0004-637X/786/2/101
  * Piaulet et al. (2021) Piaulet, C., Benneke, B., Rubenzahl, R. A., et al. 2021, AJ, 161, 70, doi: 10.3847/1538-3881/abcd3c
  * Pollacco et al. (2006) Pollacco, D. L., Skillen, I., Collier Cameron, A., et al. 2006, PASP, 118, 1407, doi: 10.1086/508556
  * Radzom et al. (2024) Radzom, B. T., Dong, J., Rice, M., et al. 2024, arXiv e-prints, arXiv:2404.06504, doi: 10.48550/arXiv.2404.06504
  * Ricker et al. (2015) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2015, Journal of Astronomical Telescopes, Instruments, and Systems, 1, 014003, doi: 10.1117/1.JATIS.1.1.014003
  * Rubenzahl et al. (2021) Rubenzahl, R. A., Dai, F., Howard, A. W., et al. 2021, AJ, 161, 119, doi: 10.3847/1538-3881/abd177
  * Sanchis-Ojeda & Winn (2011) Sanchis-Ojeda, R., & Winn, J. N. 2011, ApJ, 743, 61, doi: 10.1088/0004-637X/743/1/61
  * Siegel et al. (2023) Siegel, J. C., Winn, J. N., & Albrecht, S. H. 2023, ApJ, 950, L2, doi: 10.3847/2041-8213/acd62f
  * Skrutskie et al. (2006) Skrutskie, M. F., Cutri, R. M., Stiening, R., et al. 2006, AJ, 131, 1163, doi: 10.1086/498708
  * Sosnowska et al. (2015) Sosnowska, D., Lovis, C., Figueira, P., et al. 2015, in Astronomical Society of the Pacific Conference Series, Vol. 495, Astronomical Data Analysis Software an Systems XXIV (ADASS XXIV), ed. A. R. Taylor & E. Rosolowsky, 285, doi: 10.48550/arXiv.1509.05584
  * Southworth (2011) Southworth, J. 2011, MNRAS, 417, 2166, doi: 10.1111/j.1365-2966.2011.19399.x
  * Speagle (2020) Speagle, J. S. 2020, MNRAS, 493, 3132, doi: 10.1093/mnras/staa278
  * Stefansson et al. (2020) Stefansson, G., Mahadevan, S., Maney, M., et al. 2020, AJ, 160, 192, doi: 10.3847/1538-3881/abb13a
  * Stefansson et al. (2022) Stefansson, G., Mahadevan, S., Petrovich, C., et al. 2022, ApJ, 931, L15, doi: 10.3847/2041-8213/ac6e3c
  * Tayar et al. (2022) Tayar, J., Claytor, Z. R., Huber, D., & van Saders, J. 2022, ApJ, 927, 31, doi: 10.3847/1538-4357/ac4bbc
  * Valenti & Fischer (2005) Valenti, J. A., & Fischer, D. A. 2005, ApJS, 159, 141, doi: 10.1086/430500
  * von Zeipel (1910) von Zeipel, H. 1910, Astronomische Nachrichten, 183, 345, doi: 10.1002/asna.19091832202
  * Yee et al. (2018) Yee, S. W., Petigura, E. A., Fulton, B. J., et al. 2018, AJ, 155, 255, doi: 10.3847/1538-3881/aabfec
  * Yu & Dai (2024) Yu, H., & Dai, F. 2024, arXiv e-prints, arXiv:2406.00187, doi: 10.48550/arXiv.2406.00187
  * Zanazzi & Lai (2018) Zanazzi, J. J., & Lai, D. 2018, MNRAS, 478, 835, doi: 10.1093/mnras/sty1075