# Ponderings on the Possible Preponderance of Perpendicular Planets

Jared C. Siegel NSF Graduate Research Fellow  Department of Astrophysical
Sciences, Princeton University, 4 Ivy Lane, Princeton, NJ 08544, USA  Joshua
N. Winn Department of Astrophysical Sciences, Princeton University, 4 Ivy
Lane, Princeton, NJ 08544, USA  Simon H. Albrecht Stellar Astrophysics Centre,
Department of Physics and Astronomy, Aarhus University, Ny Munkegade 120, 8000
Aarhus C, Denmark  Department of Astrophysical Sciences, Princeton University,
4 Ivy Lane, Princeton, NJ 08544, USA

###### Abstract

Misalignments between planetary orbits and the equatorial planes of their host
stars are clues about the formation and evolution of planetary systems.
Earlier work found evidence for a peak near $90^{\circ}$ in the distribution
of stellar obliquities, based on frequentist tests. We performed hierarchical
Bayesian inference on a sample of 174 planets for which either the full three-
dimensional stellar obliquity has been measured (72 planets) or for which only
the sky-projected stellar obliquity has been measured (102 planets). We
investigated whether the obliquities are best described by a Rayleigh
distribution, or by a mixture of a Rayleigh distribution representing well-
aligned systems and a different distribution representing misaligned systems.
The mixture models are strongly favored over the single-component
distribution. For the misaligned component, we tried an isotropic distribution
and a distribution peaked at 90∘, and found the evidence to be essentially the
same for both models. Thus, our Bayesian inference engine did not find strong
evidence favoring a “perpendicular peak,” unlike the frequentist tests. We
also investigated selection biases that affect the inferred obliquity
distribution, such as the bias of the gravity-darkening method against
obliquities near $0^{\circ}$ or $180^{\circ}$. Further progress in
characterizing the obliquity distribution will probably require the
construction of a more homogeneous and complete sample of measurements.

##  1 Introduction

A planet’s orbital motion is not necessarily aligned with its host star’s
rotation. The degree of misalignment is quantified by the obliquity, $\psi$,
defined as the angle between the direction $\hat{n}_{\star}$ of the star’s
rotational angular momentum and the direction $\hat{n}_{\mathrm{o}}$ of the
planet’s orbital angular momentum. Measurements of the obliquity or its sky-
projection, $\lambda$, have been performed for several hundred transiting
exoplanets (see Triaud, 2018; Albrecht et al., 2022, for reviews).

It is hoped that such measurements will help us understand the processes of
planet formation and post-formation dynamical interactions. For example, the
observation that hot Jupiters around low-mass stars tend to have low
obliquities has been interpreted as the outcome of tidal dissipation within
the star (Winn et al., 2010), while observations of high obliquities have been
variously attributed to planet-planet scattering, Von Zeipel–Kozai–Lidov
cycles, and other mechanisms (see Section 4 of Albrecht et al., 2022). In a
few cases, there are enough clues available to implicate a particular
obliquity-excitation mechanism; for example, in the K2-290 system, there is
evidence that a companion star was responsible for tilting the protoplanetary
disk (Hjorth et al., 2021). Drawing population level conclusions will require
knowledge of the statistical distribution of stellar obliquities and its
dependence on the characteristics of planetary systems and their host stars.

Modeling of the obliquity distribution has been undertaken by Fabrycky & Winn
(2009), Morton & Winn (2014), Mazeh et al. (2015), Winn et al. (2017), Muñoz &
Perets (2018), and Louden et al. (2021), using different techniques and
progressively larger samples. In general, investigators have modeled the
distribution as a mixture of a well-aligned component and a very broad
component, and have shown evidence that the broad component consists
preferentially of “hot” stars (with masses between 1.2 and 1.9 $M_{\odot}$ and
effective temperatures between 6200 and 8000 K).

Importantly, most of the measurement techniques are not capable of determining
$\psi$ directly, but are instead sensitive to either the sky-projected
obliquity ($\lambda$) or the difference between the line-of-sight inclinations
of the spin and orbital axes ($i_{\star}$ and $i_{\rm o}$). The various angles
are related via the equation

$\hat{n}_{\star}\cdot\hat{n}_{\mathrm{o}}=\cos\psi=\cos i_{\star}\cos
i_{\mathrm{o}}+\sin i_{\star}\sin i_{\mathrm{o}}\cos\lambda.$ (1)

Recently, Albrecht et al. (2021) assembled a sample of systems for which
constraints on $\lambda$, $i_{\rm o}$, and $i_{\star}$ are all available,
usually because more than one measurement technique was applied to the same
system. In this sample, the misaligned component of the obliquity distribution
appeared to have a broad peak at 80–125∘, rather than falling off smoothly
from $0^{\circ}$ or extending all the way to 180∘. Frequentist tests assigned
low probabilities ($p\lesssim 10^{-3}$) to the hypothesis that the peak is a
statistical fluke, given the null hypothesis of isotropy.

The motivation for the work described in this paper was to (i) repeat the
frequentist tests now that more data points have been obtained, (ii) compare
different models for the obliquity distribution using Bayesian hierarchical
modeling, and (iii) investigate selection effects that might influence the
results. This paper is organized as follows. Section 2 describes the available
sample. Section 3 revisits the frequentist tests of Albrecht et al. (2021),
and Section 4 outlines the Bayesian modeling framework. The results are
presented in Section 5. In addition to analyzing the entire sample, we
examined subsamples that are expected to be influenced by tides, or for which
tides should be negligible, as described in Section 6. Selection biases are
investigated in Section 7. Our conclusions are summarized in Section 9.

##  2 Sample

Figure 1: Planetary and stellar parameters for systems with $\lambda$
measurements and the subset of systems with $\psi$ measurements (right). In
the top and bottom rows, the color of the data points conveys the planet’s
orbital separation and mass, respectively. Triangular data points are members
of the $\lambda$ sample, i.e., $\lambda$ has been measured but not $\psi$.
Circular data points are for the $\psi$ sample, i.e., both $\lambda$ and
$\psi$ have been measured.

Albrecht et al. (2022) compiled a catalog of planets for which measurements of
the sky-projected obliquity $\lambda$ are available, a subset of which also
have measurements of the full-obliquity $\psi$. We supplemented this catalog
with 16 newly available measurements (Table 1), of which 12 are for systems
that previously lacked any obliquity measurements, and 4 are for systems that
now have full-obliquity information where only the projected-obliquity
measurement was previously available. The combined catalog consists of 174
planets. Of these, 72 planets have measurements of both $\lambda$ and
$i_{\star}$; we will refer to these 72 planets as the $\psi$ sample. For the
other 102 planets, $\lambda$ has been measured but there are no available
constraints on $i_{\star}$; we will refer to these planets as the $\lambda$
sample. Figure 1 displays some key properties of the $\psi$ and $\lambda$
samples.

Most of the $\lambda$ measurements were obtained by observing the Rossiter-
McLaughlin effect, and most of the $i_{\star}$ information comes from the
combination of measurements of projected rotation velocity $v\sin i$, rotation
period, and stellar radius (Masuda & Winn, 2020). For KELT-9, KELT-17,
Kepler-13, MASCARA-4, and WASP-189 the obliquity was inferred from the
gravity-darkening method. For HAT-P-7 and HD 89345, the stellar inclination
was constrained via asteroseismology.

Table 1: Extension of the Albrecht et al. (2022) sample of projected-obliquity
$\lambda$ and obliquity $\psi$ measurements.

Name |  $\lambda$ [deg] |  $\psi$ [deg] | Update | Ref.
---|---|---|---|---
55 Cnc e | $10.0_{-20.0}^{+17.0}$ | $23.0_{-12.0}^{+14.0}$ | — | 1
HAT-P-3 b | $21.2\pm 8.7$ | $75.7_{-7.9}^{+8.5}$ | ✓ | 2
HAT-P-33 b | $5.9\pm 4.1$ | — | — | 2
HAT-P-49 b | $97.7\pm 1.8$ | — | — | 2
HD 3167 b | $6.6_{-7.9}^{+6.6}$ | $29.5_{-9.4}^{+7.2}$ | ✓ | 3
HD 3167 c | $108.9_{-5.5}^{+5.4}$ | $107.7_{-4.9}^{+5.1}$ | ✓ | 3
HD 89345 b | $74.2_{-32.5}^{+33.6}$ | $80.1_{-23.1}^{+22.3}$ | — | 2
HIP 41378 d | $57.1_{-17.9}^{+26.4}$ | — | — | 4
K2-105 b | $81.0_{-47.0}^{+50.0}$ | — | — | 2
KELT-11 b | $77.86_{-2.26}^{+2.36}$ | — | — | 5
Qatar-6 A b | $0.1\pm 2.6$ | $21.82_{-18.36}^{+8.86}$ | — | 6
TOI-1478 b | $6.2_{-5.5}^{+5.9}$ | — | — | 7
TOI-2076 b | $3.0_{-15.0}^{+16.0}$ | $18.0_{-9.0}^{+10.0}$ | — | 8
TOI-640 b | $176.0\pm 3.0$ | $104.0\pm 2.0$ | — | 9
WASP-156 b | $105.7_{-14.4}^{+14.0}$ | — | — | 2
WASP-47 b | $0.0\pm 24.0$ | $29.2_{-13.3}^{+11.1}$ | ✓ | 2

Note. — A check mark in the “Update” column indicates a planet that appeared
in the Albrecht et al. (2022) catalog for which updated measurements are now
available. References: 1 (Zhao et al., 2022), 2 (Bourrier et al., 2023), 3
(Bourrier et al., 2021), 4 (Grouffal et al., 2022), 5 (Mounzer et al., 2022),
6 (Rice et al., 2023), 7 (Rice et al., 2022), 8 (Frazier et al., 2022), 9
(Knudstrup et al., 2023).

##  3 Frequentist Hypothesis Testing

Albrecht et al. (2021) noted an apparent concentration of obliquities near
$90^{\circ}$ in the $\psi$ sample that was available at the time. Focusing
attention on the 19 planets with $\cos\psi<0.75$ in their sample, they asked:
under the null hypothesis that the spin and orbital axes are uncorrelated, how
often would the distribution of $\psi$ be concentrated near $90^{\circ}$ as
strongly as seen in the real data? The concentration was quantified by either
the dispersion of $\cos\psi$ around $0$ or by the standard deviation of
$\cos\psi$, and in both cases the null hypothesis was rejected with $p\lesssim
10^{-3}$.

In the updated sample, 25 systems satisfy the misalignment criterion
$\cos\psi<0.75$. In each of $10^{6}$ trials, we drew 25 values of $\cos\psi$
from the distribution $\mathcal{U}(-1,0.75)$ representing the null hypothesis.
We then drew values of $\cos\psi$ from the measurement posteriors of the 25
systems for which $\cos\psi<0.75$. Unlike Albrecht et al. (2021), who
approximated the measurement posteriors as skew-normal distributions, we drew
directly from the measurement posteriors. Using this method and the same
sample as Albrecht et al. (2021), the null hypothesis was rejected with median
p–values of $p=2.5\times 10^{-3}$ and $0.8\times 10^{-3}$, in terms of the
dispersion and standard deviation respectively. After adding the new and
updated measurements, the null hypothesis was rejected with median p–values of
$p=0.5\times 10^{-3}$ and $0.2\times 10^{-3}$. Thus, regardless of whether the
concentration was quantified by the dispersion or standard deviation, the null
hypothesis sample was rejected with high confidence, and the new data reduced
the corresponding $p$-values by factors of 5 and 4.

As acknowledged by Albrecht et al. (2021), these frequentist tests are
difficult to interpret because the tests were devised after seeing the data.
Also, frequentist tests can demonstrate that the data are inconsistent with
the null hypothesis, but they leave one wondering which probability
distributions are consistent with the data. Answering the latter question was
the goal of the hierarchical Bayesian modeling described below.

##  4 Hierarchical Bayesian Inference

###  4.1 Formalism

For a given dataset $\vec{d}$, a property of the system $x$ can be inferred
from the posterior probability distribution

$p(x\mid\vec{d})\propto p(\vec{d}\mid x)~{}\pi_{0}(x),$ (2)

where $p(\vec{d}\mid x)$ is the likelihood function and $\pi_{0}(x)$ is the
prior on $x$. For example, $\vec{d}$ might represent a time series of radial
velocities exhibiting the Rossiter-McLaughlin effect, and $x$ might represent
$\lambda$.

In a hierarchical model, we use a collection of $N$ datasets
$\\{\vec{d}_{1},\dots,\vec{d}_{N}\\}$ to infer the parameters $\Theta$ of a
model for the distribution $f_{x}(x\mid\Theta)$. Assuming the datasets
$\\{\vec{d}_{i}\\}$ are independent of each other, the likelihood function for
the model parameters is

$\displaystyle\mathcal{L}(\\{\vec{d}_{i}\\}\mid\Theta)$
$\displaystyle=\prod_{i=1}^{N}\int p(\vec{d}_{i}\mid
x)f_{x}(x\mid\Theta)\,dx,$ (3) $\displaystyle\propto\prod_{i=1}^{N}\int
p(x\mid\vec{d}_{i})\frac{f_{x}(x\mid\Theta)}{\pi_{0}(x)}dx.$ (4)

The population distribution $f_{x}(x\mid\Theta)$ is essentially treated as a
prior on $x$; marginalizing over $x$ then leaves a likelihood function for
$\Theta$. Throughout this work, the likelihood function for the model
parameters is denoted $\mathcal{L}(\mathrm{data}\mid\mathrm{model})$; all
other probability density functions are denoted $p(X\mid Y)$.

For our application, we modeled the obliquity distribution
$f_{\cos\psi}(\cos\psi\mid\Theta)$ using different samples of planets, and
different likelihood functions for each sample. The choices for the likelihood
function were

  * •

The $\lambda$ likelihood, based only on measurements of $\lambda$, without any
constraints on $i_{\star}$.

  * •

The $\psi$ likelihood, based on measurements of $\psi$, or separate
measurements of both $\lambda$ and $i_{\star}$.

  * •

The $\psi\,\&\,\lambda$ likelihood, i.e., the combination of $\psi$
likelihoods for systems where constraints on $\psi$ are available, and
$\lambda$ likelihoods for systems for which only $\lambda$ has been measured.

For the $\lambda$ likelihood, the data set consists of $N_{\lambda}$
posteriors for $\lambda$, and we relate $\lambda$ and $\psi$ using

$\tan\lambda=\tan\psi\sin\phi,$ (5)

an equation valid for $i_{\rm o}\approx 90^{\circ}$, as is the case for
transiting planets considered throughout this work. Here, $\phi$ is the
azimuthal angle of $\hat{n}_{\star}$ in a spherical polar coordinate system in
which $\hat{n}_{\rm o}$ is the polar axis. The azimuthal angle is assumed to
be uniformly distributed, i.e.,

$f_{\sin\phi}(x)=\frac{2}{\pi}\frac{1}{\sqrt{1-x^{2}}}.$ (6)

The likelihood function for $\Theta$ is

$\displaystyle\mathcal{L}_{\lambda}(\\{\vec{d}_{i}\\}\mid\Theta)\propto\prod_{i=1}^{N_{\lambda}}$
$\displaystyle\int_{\cos\psi,\sin\phi}d\cos\psi\times d\sin\phi$
$\displaystyle\times p_{\lambda}(\tan^{-1}[\tan\psi\sin\phi]\mid\vec{d}_{i})$
$\displaystyle\times f_{\cos\psi}(\cos\psi\mid\Theta)$ $\displaystyle\times
f_{\sin\phi}(\sin\phi),$ (7)

where $p_{\lambda}$ is the likelihood function for $\lambda$ given
$\vec{d}_{i}$, the data for the $i$th system. This likelihood function is made
computationally tractable via the $K$-samples approximation (Hogg et al.,
2010):

$\displaystyle\mathcal{L}_{\lambda}($
$\displaystyle\\{\vec{d}_{i}\\}\mid\Theta)\propto$
$\displaystyle\prod_{i=1}^{N_{\lambda}}\frac{1}{K}\sum_{k=1}^{K}\mathcal{N}(\tan^{-1}[\tan\psi_{k,\Theta}\sin\phi_{k}]\mid\lambda_{i},\sigma_{\lambda_{i}}),$
(8)

where $\cos\psi_{k,\Theta}$ is the $k$th draw from
$f_{\cos\psi}(\cos\psi\mid\Theta)$, $\sin\phi_{k}$ is the $k$th draw from
$f_{\sin\phi}$, and $K=1000$. For inference of the underlying obliquity
distribution, we found it convenient to sample from the distribution
$f_{\cos\psi}$ and approximate each measurement posterior as a normal
distribution.

For the $\psi$ likelihood, the details depend on the measurement technique. In
most cases, there are measurement posteriors for $\lambda$, $v\sin i$,
$P_{\mathrm{rot}}$, and $R_{\star}$, and the stellar inclination is
constrained using the equation

$\sin i=\frac{P_{\mathrm{rot}}}{2\pi R}~{}v\sin i,$ (9)

from which the obliquity can be calculated using Equation 1. The likelihood
function in this case is

$\displaystyle\mathcal{L}_{\psi}(\\{\vec{d}_{i}\\}\mid\Theta)\propto\prod_{i=1}^{N_{\psi}}$
$\displaystyle\int_{\cos\psi,\sin\phi,v}d\cos\psi\times d\sin\phi\times dv$
$\displaystyle\times p_{v\sin i}(v\cos\psi/\cos\lambda\mid\vec{d}_{i})$
$\displaystyle\times p_{\lambda}(\tan^{-1}[\tan\psi\sin\phi]\mid\vec{d}_{i})$
$\displaystyle\times p_{v}(v\mid\vec{d}_{i})$ $\displaystyle\times
f_{\cos\psi}(\cos\psi\mid\Theta)$ $\displaystyle\times
f_{\sin\phi}(\sin\phi).$ (10)

For notational convenience, we combined the $P_{\mathrm{rot},i}$ and
$R_{\star,i}$ measurements into a likelihood for the rotation velocity
$p_{v}(v\mid\vec{d}_{i})$. We evaluated the likelihood function via the
$K$-samples approximation, drawing from the $f_{\cos\psi}$, $f_{\sin\phi}$,
and $p_{v}$ distributions.

In the cases for which the obliquity determination was based on
asteroseismology or gravity darkening, we used the measurement posterior for
$\cos\psi$, giving a likelihood function

$\displaystyle\mathcal{L}_{\psi}(\\{\vec{d}_{i}\\}\mid\Theta)\propto\prod_{i=1}^{N_{\psi}}$
$\displaystyle\int_{\cos\psi}d\cos\psi$ $\displaystyle\times
p_{\cos\psi}(\cos\psi\mid\vec{d}_{i})$ $\displaystyle\times
f_{\cos\psi}(\cos\psi\mid\Theta).$ (11)

For our $K$-samples approximation, we approximated the posterior for
$\cos\psi$ as a normal distribution.

###  4.2 Models for the obliquity distribution

Four models for the obliquity distribution were considered. The models were
designed to determine whether the obliquity distribution has multiple
components, and if so, whether the misaligned component is peaked at a
particular value. In order of increasing complexity, the models are:

  1. 0.

A Rayleigh distribution for $\psi$ with a single parameter $\sigma_{R}$
specifying the mode and width of the distribution.

  2. 1.

A mixture of a Rayleigh distribution for $\psi$ and an isotropic distribution
(uniform in $\cos\psi$). In addition to $\sigma_{R}$, this model has a
parameter $q$ specifying the fraction of systems that are drawn from the
Rayleigh distribution.

  3. 2.

A mixture of a Rayleigh distribution for $\psi$ and an “isonormal”
distribution, a function of our devising that interpolates between an
isotropic distribution and a normal distribution (see Figure 2):

$\displaystyle f_{\mathrm{isonormal}}(\cos\psi\mid\mu,w)\propto$
$\displaystyle~{}wf_{\mathrm{isotropic}}(\cos\psi)$ $\displaystyle+$
$\displaystyle(1-w)\,\mathcal{N}\left(\cos\psi\mid\mu,\frac{w}{3}\right),$
(12)

where $\mathcal{N}$ is a normal distribution with mean $\mu$ and standard
deviation $\nicefrac{{w}}{{3}}$. As $w\rightarrow 0$, $f_{\mathrm{IG}}$ tends
toward a Dirac delta function centered at $\mu$, and as $w\rightarrow 1$,
$f_{\mathrm{isonormal}}\rightarrow f_{\mathrm{isotropic}}$. For this model,
$\mu$ is held fixed at $0$ (the “perpendicular preponderance”) and the
adjustable parameters are $\sigma_{R}$, $q$, and $w$.

  4. 3.

The same as Model 3 but allowing the mean $\mu$ of the peaked distribution to
be an adjustable parameter, for a total of four parameters:
$\\{\sigma_{R},q,\mu,w\\}$.

These models were selected because they have a minimal number of free
parameters and are easy to interpret. Each model is normalized such that
$\int_{-1}^{1}f_{\cos\psi}d\cos\psi=1$. We adopted the following priors on the
model parameters:

$\displaystyle\sigma_{R}$
$\displaystyle\sim\begin{cases}\mathcal{U}(0.3,2.0~{}\mathrm{rad})~{}\mathrm{in~{}Model~{}0},\\\
\mathcal{U}(\frac{\pi}{180},0.5~{}\mathrm{rad})~{}\mathrm{in~{}Models~{}1\mbox{--}3},\end{cases}$
$\displaystyle q$ $\displaystyle\sim\mathcal{U}(0,1),$ $\displaystyle w$
$\displaystyle\sim\mathcal{U}(0,1),$ $\displaystyle\mu$
$\displaystyle\sim\mathcal{U}(-0.5,0.5).$

The lower bounds on $\sigma_{R}$ were introduced to avoid numerical
instabilities; as long as they are sufficiently small, the specific values do
not materially affect the results.

Figure 2:  The isonormal distribution used in Model 2, for different choices
of $w$.  Table 2: Bayes Factors $R$ relative to Model 1

Restriction | Sample | $N$ | Model 0 | Model 2 | Model 3
---|---|---|---|---|---
None | $\psi$ | 72 | $(14\pm 8)$ | $\times$ | $10^{-20}$ | $2.7$ | $\pm$ | $0.5$ | $1.2$ | $\pm$ | $0.2$
| $\lambda$ | 102 | $(5\pm 2)$ | $\times$ | $10^{-14}$ | $0.382$ | $\pm$ | $0.007$ | $0.307$ | $\pm$ | $0.006$
|  $\psi$ & $\lambda$ | 174 | $(3\pm 2)$ | $\times$ | $10^{-29}$ | $1.1$ | $\pm$ | $0.8$ | $0.6$ | $\pm$ | $0.2$
No G.D. | $\psi$ | 67 | $(8\pm 6)$ | $\times$ | $10^{-20}$ | $1.6$ | $\pm$ | $0.1$ | $0.85$ | $\pm$ | $0.08$
| $\lambda$ | 102 | $(37\pm 3)$ | $\times$ | $10^{-15}$ | $0.377$ | $\pm$ | $0.008$ | $0.308$ | $\pm$ | $0.006$
|  $\psi$ & $\lambda$ | 169 | $(3\pm 4)$ | $\times$ | $10^{-30}$ | $0.5$ | $\pm$ | $0.1$ | $0.4$ | $\pm$ | $0.1$
No G.D. & Giants with Weak Tides | $\psi$ | 31 | $(34\pm 1)$ | $\times$ | $10^{-8}$ | $1.08$ | $\pm$ | $0.03$ | $0.69$ | $\pm$ | $0.02$
| $\lambda$ | 70 | $(42\pm 2)$ | $\times$ | $10^{-10}$ | $0.74$ | $\pm$ | $0.01$ | $0.434$ | $\pm$ | $0.006$
|  $\psi$ & $\lambda$ | 101 | $(76\pm 10)$ | $\times$ | $10^{-15}$ | $0.84$ | $\pm$ | $0.03$ | $0.53$ | $\pm$ | $0.02$

Note. — G.D. is gravity darkening. $R$ is the ratio of the models’ posterior
probabilities, given the data and assuming equal prior probabilities. $R>1$
implies $M_{i}$ is favored over $M_{1}$ with $R:1$ odds.

##  5 Bayesian Model Selection

The Bayes factor (or evidence ratio) was used to evaluate the competing
models. For a data set $\\{\vec{d}_{i}\\}$ and model $M$ with free-parameters
$\Theta$, the Bayesian evidence is defined as

$p(\\{\vec{d}_{i}\\}\mid M)=\int\mathcal{L}(\\{\vec{d}_{i}\\}\mid
M,\Theta)\pi_{0}(\Theta)d\Theta,$ (13)

where $\mathcal{L}(\\{\vec{d}_{i}\\}\mid M,\Theta)$ is the likelihood function
and $\pi_{0}(\Theta)$ is the prior on the model parameters. The Bayesian
evidence for each model was calculated via the $K$-samples approximation,
using a number of samples equal to the number of free parameters times
$10^{6}$. The sampling error in the calculation of the Bayesian evidence was
estimated by redrawing from the Monte Carlo samples 100 times for each model.

Given two competing models $M_{0}$ and $M_{1}$, the Bayes factor is

$R=\frac{p(\\{\vec{d}_{i}\\}\mid M_{0})}{p(\\{\vec{d}_{i}\\}\mid M_{1})}.$
(14)

If the prior probabilities of models $M_{0}$ and $M_{1}$ are equal, $R$ is the
ratio of the posterior probabilities for $M_{0}$ and $M_{1}$ given the data
$\\{\vec{d}_{i}\\}$. In this case, a Bayes factor of $R>1$ implies $M_{0}$ is
favored over $M_{1}$ with $R:1$ odds. Bayes factors of $R=10,$ $20,$ and $100$
correspond to odds ratios of $91,$ $95,$ and $99\%$, respectively; see Trotta
(2007) and Abbott et al. (2018) for reviews and applications of the Bayes
factor.

###  5.1 Full-obliquity modeling

Table 3: Parameters of the model obliquity distributions, based on complete
measurement samples.

Sample |  | $\psi$ | $\lambda$ | $\psi\,\&\,\lambda$
---|---|---|---|---
N |  | 72 | 102 | 174
Model 0 |  $\sigma_{R}$ [deg] | $44^{+3}_{-3}$ | $42^{+3}_{-3}$ | $42^{+3}_{-3}$
Model 1 |  $\sigma_{R}$ [deg] | $2.0^{+0.9}_{-0.6}$ | $9^{+3}_{-2}$ | $5^{+1}_{-1}$
| $q$ | $0.55^{+0.05}_{-0.05}$ | $0.65^{+0.06}_{-0.06}$ | $0.58^{+0.04}_{-0.03}$
Model 2 |  $\sigma_{R}$ [deg] | $2.1^{+1.0}_{-0.6}$ | $10^{+3}_{-3}$ | $5^{+1}_{-1}$
| $q$ | $0.58^{+0.05}_{-0.06}$ | $0.67^{+0.06}_{-0.06}$ | $0.6^{+0.04}_{-0.04}$
| $w$ | $0.7^{+0.2}_{-0.1}$ | $0.8^{+0.1}_{-0.2}$ | $0.8^{+0.1}_{-0.1}$
Model 3 |  $\sigma_{R}$ [deg] | $2.1^{+1.0}_{-0.6}$ | $10^{+3}_{-3}$ | $5^{+2}_{-1}$
| $q$ | $0.6^{+0.06}_{-0.06}$ | $0.66^{+0.06}_{-0.07}$ | $0.59^{+0.06}_{-0.04}$
| $w$ | $0.7^{+0.2}_{-0.2}$ | $0.8^{+0.1}_{-0.2}$ | $0.83^{+0.09}_{-0.12}$
| $\mu$ | $0.0^{+0.2}_{-0.2}$ | $0.0^{+0.3}_{-0.3}$ | $0.0^{+0.3}_{-0.2}$

Note. — Parameter posteriors for each combination of model (e.g., 0, 1, 2, or
3) and type of data ($\psi$ collection, $\lambda$ collection, and the
$\psi\&\lambda$ collection). The reported posteriors are the 16, 50, and 84th
percentiles.

Table 4: Parameters of the model obliquity distributions, based on giant
planets with weak tides and removing the gravity darkening measurements.

Sample |  | $\psi$ | $\lambda$ | $\psi\,\&\,\lambda$
---|---|---|---|---
N |  | 31 | 70 | 101
Model 0 |  $\sigma_{R}$ [deg] | $50^{+6}_{-5}$ | $40^{+3}_{-3}$ | $43^{+3}_{-3}$
Model 1 |  $\sigma_{R}$ [deg] | $2.1^{+1.7}_{-0.8}$ | $13^{+3}_{-3}$ | $8^{+2}_{-2}$
| $q$ | $0.43^{+0.08}_{-0.08}$ | $0.69^{+0.07}_{-0.08}$ | $0.56^{+0.06}_{-0.06}$
Model 2 |  $\sigma_{R}$ [deg] | $2.2^{+2.1}_{-0.9}$ | $14^{+3}_{-3}$ | $9^{+3}_{-2}$
| $q$ | $0.45^{+0.08}_{-0.08}$ | $0.73^{+0.07}_{-0.07}$ | $0.59^{+0.06}_{-0.06}$
| $w$ | $0.7^{+0.2}_{-0.2}$ | $0.7^{+0.2}_{-0.2}$ | $0.7^{+0.2}_{-0.2}$
Model 3 |  $\sigma_{R}$ [deg] | $2.1^{+2.0}_{-0.8}$ | $14^{+3}_{-3}$ | $8^{+3}_{-2}$
| $q$ | $0.45^{+0.08}_{-0.08}$ | $0.71^{+0.07}_{-0.08}$ | $0.58^{+0.06}_{-0.06}$
| $w$ | $0.7^{+0.2}_{-0.2}$ | $0.8^{+0.2}_{-0.2}$ | $0.8^{+0.1}_{-0.2}$
| $\mu$ | $0.1^{+0.2}_{-0.3}$ | $0.0^{+0.3}_{-0.2}$ | $0.1^{+0.2}_{-0.2}$

Note. — Parameter posteriors for each combination of model (e.g., 0, 1, 2, or
3) and type of data ($\psi$ collection, $\lambda$ collection, and the
$\psi\&\lambda$ collection). The reported posteriors are the 16, 50, and 84th
percentiles.

To begin, the evidence and the posteriors for the parameters of each obliquity
distribution model were inferred using the $\psi$ sample, i.e., the 72 planets
for which information about both $\lambda$ and $i_{\star}$ is available. Here
and in the other cases described below, the posteriors for the model
parameters were inferred with a Markov Chain Monte Carlo method, using 40
walkers for 50,000 iterations, of which the first third were discarded as
burn-in (a conservative threshold based on visual convergence of the walker
chains). Table 2 gives the Bayes factors relative to Model 1, and Table 3
gives the 16, 50, and 84th percentiles of the posteriors of the model
parameters.

Figure 3: Inferred obliquity distributions for different samples and
likelihoods. Regardless of the sample, the peak in the isonormal model is
small enough to be consistent or nearly consistent with the isotropic model.
Left: The $\psi$ sample of 72 planets for which measurements of both $\lambda$
and $i_{\star}$ are available, analyzed with the $\psi$ likelihoods. Center:
The $\lambda$ sample of 102 planets, analyzed with the $\lambda$ likelihoods.
Right: The entire $\psi\,\&\,\lambda$ sample of 174 planets, analyzed with the
$\psi\,\&\,\lambda$ likelihoods. Solid lines demarcate the mean, and the
shaded regions show the $1\sigma$ uncertainty intervals.

In Model 0 (Rayleigh), the width $\sigma_{R}$ was found to be $44^{\circ}\pm
3^{\circ}$. This model was strongly disfavored ($R<10^{-9}$) relative to the
other models in which the misaligned systems were described by a separate
component. We conclude that the obliquity distribution is best described as a
mixture between a well-aligned component and a misaligned component, rather
than a single smooth distribution.

In the multi-component models (1–3), the fraction of systems belonging to the
Rayleigh component was $\approx$60%, and $\sigma_{R}$ was found to be
$\approx$2∘. Regarding the misaligned components, the Bayesian evidence favors
the peaked models (2 and 3) over the isotropic distribution (1), but not
significantly. Model 2 (Rayleigh $+$ 90∘ peak) is favored over Model 1
(Rayleigh $+$ isotropic) with $R\approx 3$. Model 3, in which the peak
location is allowed to vary, is consistent with a peak at 90∘; however, the
constraint on $\mu$ is weak and the degree of concentration is low. For Models
2–3, the width parameter $w$ was found to be $\gtrsim$0.7, indicating a broad
distribution. The inferred obliquity distributions are visualized in Figure 3.

For reference, we repeated the inference using the same sample of 72 planets,
but using only the $\lambda$ likelihoods and ignoring the constraints on the
stellar inclination angle. The results are nearly the same as those that were
obtained using the $\psi$ likelihoods, all consistent within $1\sigma$.
Evidently, the inclination information has only a minor effect on the
results.111This conclusion was also reached by Dong & Foreman-Mackey in work
shared with us prior to publication.

###  5.2 Projected-obliquity modeling

Next, we considered the $\lambda$ sample, consisting of 102 planets for which
only $\lambda$ was measured and no information is available about the stellar
inclination. Table 2 gives the Bayes factors relative to Model 1, and Table 3
gives the 16, 50, and 84th percentiles of the model posteriors. The inferred
obliquity distributions are visualized in Figure 3.

Again, Model 0 (Rayleigh) was strongly disfavored relative to the multi-
component models. Also similar to the previous case is that Models 1–3 were
found to be effectively indistinguishable by the Bayesian evidence; Model 1 is
favored with odds ratios $\approx$3\. In Models 2–3, the width parameter $w$
was found to be $0.8_{-0.2}^{+0.1}$. Relative to the $\psi$ sample, the
$\lambda$ sample favors a larger width for the well-aligned component
($\sigma_{R}\approx 10^{\circ}$ compared to $\sigma_{R}\approx 2^{\circ}$) and
provides weaker evidence for models with a peaked misaligned component; these
trends may indicate underlying differences between the $\psi$ and $\lambda$
sample, which we investigate in Section 7.

###  5.3 Joint obliquity and projected-obliquity modeling

Finally, we took advantage of all the available information. We combined the
collection of $\psi$ constraints from the 72 systems in the $\psi$ sample with
the $\lambda$ constraints from the 102 systems in the $\lambda$ sample. The
Bayes factors relative to Model 1 are presented in Table 2. Table 3 presents
the 16, 50, and 84th percentiles of the model posteriors. The inferred
obliquity distributions are visualized in Figure 3. Model 0 was strongly
disfavored relative to Models 1–3, and Models 1–3 were effectively
indistinguishable. The preferred models are a mixture of a narrow Rayleigh
distribution and a very broad or isotropic component.

To summarize, the single-component Model 0 was always trounced by the multi-
component Models 1–3. While the results of analyzing the $\psi$ and $\lambda$
systems differ in detail, neither the isotropic model nor the peaked models
were strongly favored by the Bayesian evidence.

Figure 4:  Projected-obliquity (left) and full-obliquity (right) versus
stellar effective temperature. Blue points represent planets for which
$\lambda$ has been measured but not $\psi$. Purple points represent planets
for which $\psi$ has been measured. Circles are the giant-planet hosts in the
weak-tides sample, triangles are the gravity-darkening measurements, and
crosses are all other measurements.

##  6 Sample dependence

The sample of planets covers a wide range of orbital separations, planet
masses, and stellar types. The resulting diversity of possible modes of planet
formation and evolution makes the obliquity distribution of the entire sample
difficult to interpret. In particular, for massive close-orbiting planets,
dissipative tidal interactions generally lead to orbital circularization,
synchronization of spin and orbital periods, and spin-orbit alignment. Tidal
dissipation is expected to be faster for low-mass host stars with thick
convective envelopes (with effective temperatures below the “Kraft break” at
about 6250 K) than for more massive stars with radiative envelopes. Although
the mechanism by which tidal dissipation occurs is uncertain, the obliquity
data shows the qualitative trends one would expect from tidal dissipation.
Close-orbiting giant planets around low-mass stars generally have low
obliquities, while a broader distribution of obliquities is seen in systems
with smaller planets, wider orbits, or hotter stars (Winn et al., 2010;
Albrecht et al., 2012, 2022).

In an effort to simplify the interpretation of the obliquity distribution, we
repeated our Bayesian model comparisons after applying restrictions to the
samples intended to group together systems where one might expect similar
formation mechanisms and obliquity-altering mechanisms to have occurred. We
focused on giant planets ($M_{p}>0.3M_{\rm Jup})$, and included those with
$a/R_{\star}<300$ (a criterion that excluded only $\beta$ Pic b). We omitted
the 5 systems for which the obliquity was determined with the gravity
darkening method, out of concern about a selection bias favoring large
misalignments (see Section 7.4). For a system to be categorized as having
“weak tides”, we required $a/R_{\star}>7$ or $T_{\mathrm{\rm eff}}<6250$ K.
Figure 4 shows key properties of the samples and identifies the weak-tides
systems.

Figure 5:  The $\lambda$ measurements from the $\psi$ sample are slightly but
systematically higher than the $\lambda$ measurements from the $\lambda$
sample. This effect is likely due at least in part to the selection bias of
the gravity darkening method, which was used in 5 cases in the $\psi$ sample
and no cases in the $\lambda$ sample. In each panel, the $\lambda$
measurements are plotted from lowest to highest, with blue points for the
$\lambda$ sample and purple points for the $\psi$ sample. The left panel shows
all available measurements. The center panel omits data from the gravity
darkening method. The right panel is for the weak-tides sample. For each
measurement, the $1\sigma$ uncertainty interval is plotted as a vertical line.

###  6.1 Weak tides

After applying the weak-tides restriction, the $\lambda$ sample consists of 70
giant planets for which only the projected obliquity has been measured, and
the $\psi$ sample consists of 31 giant planets for which full-obliquity
information is available. Because tidal dissipation is expected to be
negligible, the obliquity distribution inferred for these samples may more
closely reflect the outcome of the formation and migration mechanisms.

As in Section 5, the parameters were inferred for 4 different models for the
obliquity distribution, and different choices for the sample and likelihood
function. Table 4 gives the credible intervals of the parameters of the models
for the obliquity distribution.

The results of analyzing the weak-tides sample were largely consistent with
the results of analyzing the unrestricted sample. We expected the weak-tides
sample to give a lower value of $q$, the fraction of systems in the well-
aligned Rayleigh component, and a larger value of $\sigma_{R}$, the scale
factor of the well-aligned component, but these differences did not emerge at
a statistically significant level. Instead, discrepancies arose between the
model posteriors inferred from the $\lambda$ and $\psi$ samples. Relative to
the $\psi$ sample, the $\lambda$ sample favored a larger width for the well-
aligned component ($\sigma_{R}\approx 10^{\circ}$ vs. $2^{\circ}$) and a
higher fraction of well-aligned planets ($q\approx 0.7$ vs. 0.5).

Table 2 gives the Bayes factor of each model (relative to Model 1) for samples
with weak tides. Regardless of the sample and likelihoods that were used, the
single Rayleigh distribution was strongly disfavored relative to the multi-
component models, and the multi-component models were supported by nearly
equivalent evidence.

###  6.2 Strong tides

We also constructed a sample of giant planets for which tides are expected to
be strong, using the criteria $M_{p}>0.3M_{\rm Jup}$, $a/R_{\star}>7$, and
$T_{\mathrm{\rm eff}}<6250$ K. After applying these restrictions, the $\psi$
sample has 11 members and the $\lambda$ sample has 12 members. These systems
might be expected to have systematically low obliquities due to tidal
obliquity damping, and indeed, the data are all consistent with zero
obliquity, with varying degrees of uncertainty.

For this sample, the multi-component models are unjustified, and only Model 0
was fitted. For the three choices of sample and likelihood ($\psi$, $\lambda$,
and $\psi\,\&\,\lambda$), the Rayleigh scale factor was found to be
$4^{+3}_{-2}$, $5^{+4}_{-2}$, and $4^{+2}_{-2}[^{\circ}]$, respectively. The
$95$th percentile upper bounds on the Rayleigh scale factor are $9$, $11$, and
$7[^{\circ}]$, respectively.

##  7 Selection biases

The frequentist tests of the $\psi$ sample assigned a low probability to an
isotropic distribution for the misaligned systems, and our Bayesian analysis
of the $\psi$ sample gave somewhat stronger evidence for a $90^{\circ}$ peak
than an isotropic distribution. Yet, our analysis of the $\lambda$ sample
favored an isotropic distribution. Additionally, the posterior distributions
of the model parameters only partially overlapped between the $\psi$ and
$\lambda$ samples. Although the Bayes factors were modest in all cases, we
wondered how and why the $\psi$ sample might differ from the $\lambda$ sample.
Figure 5 shows the projected-obliquity measurements from the $\lambda$ sample
alongside those from the $\psi$ sample. In this figure, the $\lambda$
measurements in the $\psi$ sample appear to be systematically higher than
those in the $\lambda$ sample.

This made us question some premises implicit in all of the work described
above, namely, that the planets with full-obliquity information are drawn
randomly from the same parent population as the planets with only projected-
obliquity information, and that the availability of full-obliquity information
is independent of the value of the obliquity. In this section, we consider the
possibility that these premises are false. Below, in Section 7.1, we describe
some statistical tests, and in Sections 7.2, 7.3, and 7.4, we discuss some
possible sources of bias inherent in the measurement techniques.

Figure 6:  The distribution of achievable measurement precision in $\lambda$
for the nominal system described in the text, assuming random stellar
orientations and either a fixed impact parameter or a uniform distribution of
impact parameters. The median value of $\sigma_{\lambda}$ is presented as a
solid line, and the shaded regions demarcate the $1\sigma$ and $2\sigma$
percentiles. Color indicates the assumed radial velocity measurement
uncertainty: $1~{}\mathrm{m~{}s}^{-1}$ (blue) or $5~{}\mathrm{m~{}s}^{-1}$
(gold).

###  7.1  Observational evidence

To test the null hypothesis that the stellar obliquities in the $\psi$ and
$\lambda$ samples are drawn from the same distribution, we applied the
Kolmogorov-Smirnov (KS) test in a series of 1,000 Monte Carlo trials. For a
given trial, we drew (with replacement) $N_{\psi}$ values of $\lambda$ from
the $\psi$ sample, and $N_{\lambda}$ values of $\lambda$ from the $\lambda$
sample. These values were drawn from the measurement posteriors for each
system. Then, the two-sample KS test was applied to compute the $p$–value for
the null hypothesis that the $\lambda$ values are drawn from the same
distribution. The null hypothesis was “ruled out” ($p<5\%$) in $\lesssim 5\%$
of the trials. This was true when the method was applied to the entire planet
sample, and when it was applied after enforcing either the weak-tides or
strong-tides restrictions. Two-sample Anderson-Darling tests also ruled out
the null hypothesis in $\lesssim 5\%$ of the trials, for each choice of sub-
sample. The median $p$-value was $50\%$ for the whole sample and $40\%$ for
the tidally-decoupled sample. Thus, these tests did not reject the hypothesis
that the $\lambda$ measurements in the $\psi$ sample are drawn from the same
distribution as the $\lambda$ measurements in the $\lambda$ sample.

###  7.2  Rossiter-McLaughlin effect

For most of the systems in the $\psi$ sample, $\lambda$ was measured with the
Rossiter-McLaughlin (RM) effect and $i_{\star}$ was measured or constrained
with a different method. The RM effect, in brief, is the distortion of a
star’s spectral lines that occurs during transits due to stellar rotation. For
example, when a planet blocks a portion of the approaching half of the stellar
disk, the blueshifted components of the spectral lines are slightly
suppressed, which can be detected directly in the line profile or as an
apparent redshift of the entire line. This “anomalous” Doppler shift depends
chiefly on $v\sin i$, $\lambda$, the transit depth $\delta\approx(r/R)^{2}$,
and the transit impact parameter $b$. Gaudi & Winn (2007) derived an
approximate formula, later corrected by Albrecht et al. (2022), for the
achievable uncertainty in $\lambda$ as a function of the system parameters:

$\sigma_{\lambda}=\frac{\sigma_{v}/\sqrt{N}}{v\sin
i}\left(\frac{r}{R}\right)^{-2}\left[\frac{(1-b^{2})\cos^{2}\lambda+3b^{2}\sin^{2}\lambda}{b^{2}(1-b^{2})}\right]^{1/2}.$
(15)

The derivation assumes there are $N$ uniformly spaced radial velocity
measurements, each with an uncertainty of $\sigma_{v}$. Limb darkening is
neglected, as are the uncertainties in $b$ and $r/R$.

The dependence of $\sigma_{\lambda}$ on $\lambda$ raises the question of
whether the sample of systems with published obliquity measurements is biased
in favor of obliquities that tend to allow for smaller values of
$\sigma_{\lambda}$ and consequently higher probabilities to detect the RM
effect. To understand this source of bias, we conducted a series of Monte
Carlo trials involving a hypothetical 1 $R_{\odot}$ star with rotation
velocity $v=2.5~{}\mathrm{km~{}s}^{-1}$ and a transiting Jupiter-sized planet.
In each trial, the star’s orientation was drawn randomly from an isotropic
distribution. Figure 6 shows $\sigma_{\lambda}$ as a function of $\psi$ for
three illustrative values of the impact parameter, and for two choices of
$\sigma_{v}$. The bias depends strongly on $b$, reflecting the limiting
behaviors of Eqn 15: $\sigma_{\lambda}\propto\cos^{2}\lambda$ as $b\rightarrow
0$, and $\sigma_{\lambda}\propto\sin^{2}\lambda$ as $b\rightarrow 1$. The
rightmost panel of Figure 6 shows that after marginalizing over a uniform
distribution of $b$ between 0 and 1, $\sigma_{\lambda}$ varies weakly with
$\psi$. Values of $\psi=0^{\circ}$ and $180^{\circ}$ lead to average values of
$\sigma_{\lambda}$ that are about two times lower than when $\psi=90^{\circ}$.
Thus, we expect that any such bias in the real sample to be small and to favor
well-aligned and anti-aligned systems rather than producing a peak at
$\psi=90^{\circ}$.

Figure 7:  Biases in the $v\sin i$ method (Section 7.3) and the gravity
darkening method (Section 7.4) shape the recovered obliquity distribution.
Each panel shows the obliquity distribution of hypothetical stars drawn from a
randomly-oriented population, after applying increasingly stringent detection
thresholds. Left: Obliquities measured with the $v\sin i$ method, using three
detection thresholds: (i) $v\sin i>2~{}\mathrm{km~{}s}^{-1}$, (ii) photometric
amplitude $A>0.7\%$, where $A=f_{\mathrm{spot}}\sin i$, and (iii) $A>0.9\%$.
The $v\sin i$ method favors $\psi=0^{\circ}$ and $180^{\circ}$. Right:
Obliquities measured with the gravity-darkening method. The detection
thresholds were parameterized by the difference in reduced-$\chi^{2}$ between
modeling the light curve with gravity-darkening or a standard transit profile.

###  7.3 The $v\sin i$ method

Many of the constraints on stellar inclination come from the combination of
measurements of the star’s projected rotation velocity $v\sin i$, radius $R$,
and rotation period $P_{\mathrm{\rm rot}}$, based on Equation 9. Any selection
biases regarding $R$ are unlikely to depend on the star’s orientation, but
$v\sin i$ and $P_{\rm rot}$ are potentially problematic. Measurements of
$v\sin i$ tend to be limited to an accuracy of 1–2 km s-1 by systematic
effects; thus, we would expect the fractional uncertainty in $v\sin i$ to be
larger when $\sin i$ is small. Small values of $\sin i$ also make measurements
of $P_{\rm rot}$ more difficult, because the photometric variability due
starspots and plages has a lower amplitude at low inclinations.

We conducted more Monte Carlo trials to understand the consequent selection
biases as a function of $\psi$. The intrinsic stellar properties were held
fixed: $R=R_{\odot}$, $P_{\mathrm{\rm rot}}=15$ days, and a maximum loss of
light due to spots of $f_{\rm spot}=1\%$. In each trial, the star’s
orientation and the planet’s obliquity were drawn independently from an
isotropic distribution. We assumed the photometric variability amplitude $A$
to be $f_{\rm spot}\sin i$, an approximation that assumes the spots are small
and nearly equatorial, and neglects the effects of limb darkening (see, e.g.,
Jackson & Jeffries, 2012; Mazeh et al., 2015). The stellar inclination was
considered “measurable” when $v\sin i>2~{}\mathrm{km~{}s}^{-1}$ and $A>A_{\rm
min}$, or equivalently,

$\sin
i>\max\left[\frac{2~{}\mathrm{km~{}s}^{-1}}{v},\frac{A_{\min}}{f_{\mathrm{spot}}}\right].$
(16)

For the adopted stellar parameters, the $v\sin i$ criterion is the limiting
factor when $A_{\min}\lesssim 0.6\%$. Figure 7 shows the obliquity
distribution of stars with measurable inclinations, for increasingly stringent
detection criteria.

These selection effects would cause the sample to be biased in favor of well-
aligned and anti-aligned systems, for which the inclination is guaranteed to
be close to $90^{\circ}$. Misaligned systems, and perpendicular systems in
particular, would occasionally be observed pole-on and are thereby disfavored.

###  7.4  Gravity Darkening

The photosphere of a rapidly rotating star is cooler and fainter near the
equator than the poles, because the equatorial region is centrifugally lifted
to higher elevation. This effect, known as gravity darkening, causes
distortions of a transit light curve that depend on both $\lambda$ and
$i_{\star}$ in addition to the usual transit parameters (Barnes, 2009). For
five of the systems in the $\psi$ sample (KELT-9, KELT-17, Kepler-13,
MASCARA-4, and WASP-189), the stellar obliquity was determined by modeling
transit light curves affected by gravity darkening. In contrast, none of the
measurements in the $\lambda$ sample relied on gravity darkening. Could this
difference account for the differences in the inferred obliquity distributions
of the $\psi$ and $\lambda$ samples?

Gravity-darkening anomalies in the light curve are typically on the order of
$\lesssim$0.1%, making them challenging to detect, and the detectability
depends on the stellar orientation. When $\psi$ is near $0^{\circ}$ or
$180^{\circ}$, the transit light curve retains its usual “U-shape” with mirror
symmetry around the time of minimum light, causing the parameters describing
gravity darkening to be strongly covariant with the usual transit parameters.
For other obliquities, the transit light curve develops an asymmetry or a
“W-shape”, either of which reduces the covariance between the gravity-
darkening parameters and the other parameters, allowing the gravity-darkening
signal to be detected and interpreted with greater confidence. Thus, the
sample of systems for which the gravity darkening signal has been reported is
likely to be biased against well-aligned and anti-aligned systems.

To understand the extent of this bias, we performed another round of Monte
Carlo trials involving a single star that can have any obliquity and a
transiting planet that can have any impact parameter. Based on the properties
of one of the best studied systems, Kepler-13A (Masuda, 2015), we fixed the
properties of the hypothetical star to be $M=1.80\,M_{\odot}$,
$R=1.71\,R_{\odot}$, $P_{\mathrm{\rm rot}}=20$ hours, $r/R=0.085$, and $P_{\rm
orb}=1.76$ days. A quadratic limb-darkening law was adopted, with coefficients
given by Masuda (2015): $u_{1}+u_{2}=0.53$ and $u_{1}-u_{2}=0.26$. The
effective temperature $T$ was assumed to vary with surface gravity $g$ as

$T=(7500\,{\rm K})\left(\frac{g}{g_{\mathrm{pole}}}\right)^{1/4}.$ (17)

In each of $2\times 10^{5}$ trials, the stellar orientation was drawn from an
isotropic distribution, and the planet’s transit impact parameter was drawn
from a uniform distribution between 0 and 1. We generated a transit light
curve exhibiting the effects of gravity darkening using the pytransit code
written by Parviainen (2015). We calculated the difference in
reduced-$\chi^{2}$ between the best-fit model including gravity darkening and
the best-fit model with no gravity darkening. For computational efficiency,
only $i_{\star}$, $\lambda$, $P_{\mathrm{rot}}$, the limb-darkening
coefficients, and $r/R$ were treated as free parameters. Through these Monte
Carlo trials, we mapped the relation between $\Delta\chi^{2}/\nu$ and $\psi$.

Figure 7 shows the obliquity distribution of those systems for which
$\Delta\chi^{2}/\nu$ exceeds a threshold value. Results for three different
choices of the threshold value are shown. Even for modest detection thresholds
(i.e., $\Delta\chi^{2}/\nu>2$), the obliquity distribution of the systems with
detectable gravity-darkening signals is biased strongly against
$\psi=0^{\circ}$ and $180^{\circ}$.

The bias goes in the right direction to explain why the $\psi$ sample shows
stronger evidence for a “pile-up” of obliquities near $90^{\circ}$ than the
$\lambda$ sample. To pursue this point, we repeated the modeling procedure
described in Section 5 after excluding the 5 gravity-darkening systems from
the sample. The resulting Bayes factors, relative to Model 1, are presented in
the third row of Table 2. After omitting the gravity-darkening measurements,
the $\psi$ sample shows a weaker preference for a peak at 90∘ in the obliquity
distribution. This trend is evident in Figure 5.

Figure 8:  Observational investment required to recover a peaked obliquity
distribution, for three idealized campaigns: (i) measuring $\lambda$ only,
(ii) measuring $\sin i$ only, and (iii) measuring both $\lambda$ and $\sin i$.
The vertical axis is the fraction of realizations in which the postulated peak
(see text) is successfully identified. In the left panel, the horizontal axis
is the sample size; in the right panel, the horizontal axis is the total
observing time, assuming $5$ hr per $\lambda$ measurement and $0.25$ hr per
$v\sin i$ measurement.

##  8 Population Synthesis

Bayesian analysis of the $\psi$ and $\lambda$ samples did not distinguish
between a peaked or flat distribution for the misaligned systems. While this
result held for both the $\psi$ and $\lambda$ samples—as well as the weak-
tides sub-sample—minor discrepancies persist: the $\psi$ sample marginally
favors a $90^{\circ}$ peak, more so than the $\lambda$ sample, and the
posteriors of the model parameters only partially overlap between the $\psi$
and $\lambda$ samples. Prompted by these mysteries, we investigated selection
biases in Section 7.

However, the impacts of selection biases are difficult to constrain for
heterogeneous samples. Stars with different characteristics may have different
stellar obliquity distributions, and they may also have different
characteristics that affect whether the stellar obliquity can be measured. For
example, the ability to measure rotation periods from photometry is easier for
cool spotted stars with convective envelopes than for hotter stars with
radiative envelopes, and there is also evidence that cool and hot stars have
intrinsically different obliquity distributions. We advocate for the creation
of a more homogeneous sample to understand and overcome these biases. To
achieve this goal, would it be better to focus on $\lambda$ measurements,
which are observationally expensive, or on $\sin i$ measurements, which bear
less information but are less expensive?

We performed some Monte Carlo experiments to begin answering this question. We
generated a synthetic sample of obliquities drawn from an isonormal model with
$w=0.7$, a level of “peakedness” motivated by the fits to the current sample.
Host star properties and measurement uncertainties were assigned for each draw
by randomly selecting host stars and their measurement uncertainties from the
observed sample. This process neglected the dependence of $\sigma_{\lambda}$
on $\lambda$ and $\sigma_{v\sin i}$ on $v\sin i$, but retained the $v\sin i$
selection bias by enforcing realistic uncertainties $\sigma_{v\sin i}$. For a
given synthetic sample, we computed the Bayes factor comparing the isonormal
model and an isotropic model. The peaked distribution was deemed “recovered”
if the Bayes factor for the isonormal model relative to the isotropic model
was greater than $10$. When calculating Bayes factors, we considered three
scenarios: (i) measuring $\lambda$ only, via the RM effect; (ii) measuring
$\sin i$ only, by combining $v\sin i$, $R_{\star}$, and $P_{\mathrm{rot}}$;
and (iii) measuring both $\lambda$ and $\sin i$ (the $\psi$ sample). We did
not consider gravity darkening measurements because of the strong bias of that
method against aligned and anti-aligned stars. For a given choice of sample
size and measurement technique, 500 trials were conducted.

Figure 8 shows the fraction of trials in which the peaked model was
successfully recovered. To reach a success rate of $>90\%$, the necessary
sample size was 180, 200, and 1000, for the $\psi$, $\lambda$, and $\sin i$
samples, respectively. As expected from the results presented in Section 5.1,
a single $\lambda$ measurement is typically more valuable than a single $\sin
i$ measurement. However, a typical RM measurement requires a time-critical
observation lasting $\approx\,5$ hr on a large telescope with a stabilized
high-resolution spectrograph. In contrast, measuring $v\sin i$ with such a
telescope requires no more than 15 minutes (including overheads), does not
require a highly stable wavelength calibration, and can be conducted at any
time. Applying these fiducial costs of 5 hr per $\lambda$ measurement and
$0.25$ hr per $\sin i$ measurement, we found that the necessary observing time
is 990, 980, and 250 hr for the $\psi$, $\lambda$, and $\sin i$ samples,
respectively. Thus, in these simulations, the $v\sin i$ technique is favored,
even without accounting for bad weather and scheduling constraints that would
also favor the $v\sin i$ technique. Because of all the simplifying
assumptions, these results should be considered illustrative rather than
quantitatively reliable.

##  9 Summary and Conclusions

Even before the first misaligned planetary systems were discovered, theorists
explored mechanisms for tilting a planet’s orbital plane away from the star’s
equatorial plane (e.g., Yu & Tremaine, 2001; Fabrycky & Tremaine, 2007).
Theoretical activity accelerated after the first few sightings of grossly
misaligned systems (e.g., Hébrard et al., 2008; Winn et al., 2009; Queloz et
al., 2010), leaving us with numerous proposed mechanisms for spin-orbit
misalignment, such as magnetic warping of the protoplanetary disk (Foucart &
Lai, 2011), Von Zeipel–Kozai–Lidov cycles (Fabrycky & Tremaine, 2007; Naoz et
al., 2011), planet-planet scattering (Chatterjee et al., 2008), and
inclination excitation from secular resonance crossing (Petrovich et al.,
2020). Characterizing the observed obliquity distribution is an important goal
in testing these theories for spin-orbit misalignment.

We applied hierarchical Bayesian inference to a sample of 174 planets (72 with
full-obliquity constraints and 102 with only projected-obliquity constraints).
Through Bayesian model selection, we tried to establish: (i) Is the obliquity
distribution compatible with a single smooth distribution centered on
$0^{\circ}$? (ii) Are the misaligned systems consistent with isotropy or do
they exhibit a concentration? (iii) If so, what is the preferred value of the
obliquity?

A single Rayleigh distribution was found to be a terrible description of the
data, relative to the models in which some systems are drawn from a Rayleigh
distribution and others are drawn from a broader distribution. The sample with
full-obliquity constraints showed greater evidence for a peak in the obliquity
distribution at about $90^{\circ}$, relative to a model in which the
misaligned systems are drawn from an isotropic distribution. However, the odds
ratios were modest ($\lesssim$3), and the results are likely affected by the
selection bias of the gravity-darkening method that favors near-perpendicular
systems. Analyzing the larger sample of systems with only projected-obliquity
constraints, or the union of all systems with any kind of obliquity
constraints, the Bayesian evidence did not reveal any preference for a
concentration near 90∘. Our investigation of selection biases showed that the
gravity-darkening method favors systems with large obliquities, as noted
above, while the other methods favor aligned and anti-aligned systems.

We believe that major progress on measuring the obliquity distribution will
require the construction of a large sample of systems that have been measured
with the same technique (or combination of techniques), coupled with a good
understanding of the selection function. It is not clear how to achieve this
goal. The RM technique is observationally expensive, typically requiring 4-5
hours per system using a large telescope equipped with a high-resolution
spectrograph. The gravity-darkening technique requires transiting planets
around rapidly rotating stars, of which few are known; it also requires very
precise transit photometry. The $v\sin i$ technique is the only one for which
it seems feasible to expand the sample to thousands of stars or more. Only a
single high-resolution spectrum of each system is needed, along with the
stellar rotation period, which can be obtained from ongoing or planned wide-
field photometric time-domain surveys. The ideal types of star are probably
early G stars or late F stars, for which moderately rapid rotation is expected
($\gtrsim$ 5 km/s) and for which rotationally-induced photometric variations
are also expected. Our population-synthesis calculations demonstrated that a
$v\sin i$ survey would require less observing time than a RM survey (by a
factor of a few) to distinguish between a peaked and isotropic obliquity
distribution, although these calculations were based on many simplifying
assumptions. Before deciding on an observing program, other considerations
would need to be taken into account, such as the availability and reliability
of photometric rotation periods, and the limited number of transiting planets
that are known to exist around early G and late F stars. Nobody said that
solving the spin-orbit misalignment problem would be easy.

We thank Jiayin Dong and Dan Foreman-Mackey for helpful conversations and for
sharing their preprint of a paper about modeling the obliquity distribution.
JS acknowledges support by the National Science Foundation Graduate Research
Fellowship Program under Grant DGE-2039656. SHA acknowledges the hospitality
of the Department of Astrophysical Sciences at Princeton University and
support from the Danish Council for Independent Research through a grant, No.
2032-00230B. Any opinions, findings, and conclusions or recommendations
expressed in this material are those of the author(s) and do not necessarily
reflect the views of the National Science Foundation.

## References

  * Abbott et al. (2018) Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2018, Phys. Rev. D, 98, 043526, doi: 10.1103/PhysRevD.98.043526
  * 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
  * Barnes (2009) Barnes, J. W. 2009, ApJ, 705, 683, doi: 10.1088/0004-637X/705/1/683
  * 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. (2023) Bourrier, V., Attia, O., Mallonn, M., et al. 2023, A&A, 669, A63, doi: 10.1051/0004-6361/202245004
  * Chatterjee et al. (2008) Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580, doi: 10.1086/590227
  * Fabrycky & Tremaine (2007) Fabrycky, D., & Tremaine, S. 2007, ApJ, 669, 1298, doi: 10.1086/521702
  * Fabrycky & Winn (2009) Fabrycky, D. C., & Winn, J. N. 2009, ApJ, 696, 1230, doi: 10.1088/0004-637X/696/2/1230
  * Foucart & Lai (2011) Foucart, F., & Lai, D. 2011, MNRAS, 412, 2799, doi: 10.1111/j.1365-2966.2010.18176.x
  * Frazier et al. (2022) Frazier, R. C., Stefansson, G., Mahadevan, S., et al. 2022, arXiv e-prints, arXiv:2212.06266, doi: 10.48550/arXiv.2212.06266
  * Gaudi & Winn (2007) Gaudi, B. S., & Winn, J. N. 2007, ApJ, 655, 550, doi: 10.1086/509910
  * Grouffal et al. (2022) Grouffal, S., Santerne, A., Bourrier, V., et al. 2022, A&A, 668, A172, doi: 10.1051/0004-6361/202244182
  * Hébrard et al. (2008) Hébrard, G., Bouchy, F., Pont, F., et al. 2008, A&A, 488, 763, doi: 10.1051/0004-6361:200810056
  * Hjorth et al. (2021) Hjorth, M., Albrecht, S., Hirano, T., et al. 2021, Proceedings of the National Academy of Science, 118, e2017418118, doi: 10.1073/pnas.2017418118
  * Hogg et al. (2010) Hogg, D. W., Myers, A. D., & Bovy, J. 2010, ApJ, 725, 2166, doi: 10.1088/0004-637X/725/2/2166
  * Jackson & Jeffries (2012) Jackson, R. J., & Jeffries, R. D. 2012, MNRAS, 423, 2966, doi: 10.1111/j.1365-2966.2012.21119.x
  * Knudstrup et al. (2023) Knudstrup, E., Albrecht, S. H., Gandolfi, D., et al. 2023, arXiv e-prints, arXiv:2302.01702, doi: 10.48550/arXiv.2302.01702
  * Louden et al. (2021) Louden, E. M., Winn, J. N., Petigura, E. A., et al. 2021, AJ, 161, 68, doi: 10.3847/1538-3881/abcebd
  * Masuda (2015) Masuda, K. 2015, ApJ, 805, 28, doi: 10.1088/0004-637X/805/1/28
  * Masuda & Winn (2020) Masuda, K., & Winn, J. N. 2020, AJ, 159, 81, doi: 10.3847/1538-3881/ab65be
  * Mazeh et al. (2015) Mazeh, T., Perets, H. B., McQuillan, A., & Goldstein, E. S. 2015, ApJ, 801, 3, doi: 10.1088/0004-637X/801/1/3
  * Morton & Winn (2014) Morton, T. D., & Winn, J. N. 2014, ApJ, 796, 47, doi: 10.1088/0004-637X/796/1/47
  * Mounzer et al. (2022) Mounzer, D., Lovis, C., Seidel, J. V., et al. 2022, A&A, 668, A1, doi: 10.1051/0004-6361/202243998
  * Muñoz & Perets (2018) Muñoz, D. J., & Perets, H. B. 2018, AJ, 156, 253, doi: 10.3847/1538-3881/aae7d0
  * Naoz et al. (2011) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2011, Nature, 473, 187, doi: 10.1038/nature10076
  * Parviainen (2015) Parviainen, H. 2015, MNRAS, 450, 3233, doi: 10.1093/mnras/stv894
  * 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
  * Queloz et al. (2010) Queloz, D., Anderson, D. R., Collier Cameron, A., et al. 2010, A&A, 517, L1, doi: 10.1051/0004-6361/201014768
  * Rice et al. (2023) Rice, M., Wang, S., Gerbig, K., et al. 2023, AJ, 165, 65, doi: 10.3847/1538-3881/aca88e
  * Rice et al. (2022) Rice, M., Wang, S., Wang, X.-Y., et al. 2022, AJ, 164, 104, doi: 10.3847/1538-3881/ac8153
  * Triaud (2018) Triaud, A. H. M. J. 2018, in Handbook of Exoplanets, ed. H. J. Deeg & J. A. Belmonte, 2, doi: 10.1007/978-3-319-55333-7_2
  * Trotta (2007) Trotta, R. 2007, MNRAS, 378, 72, doi: 10.1111/j.1365-2966.2007.11738.x
  * Winn et al. (2010) Winn, J. N., Fabrycky, D., Albrecht, S., & Johnson, J. A. 2010, ApJ, 718, L145, doi: 10.1088/2041-8205/718/2/L145
  * Winn et al. (2009) Winn, J. N., Johnson, J. A., Fabrycky, D., et al. 2009, ApJ, 700, 302, doi: 10.1088/0004-637X/700/1/302
  * Winn et al. (2017) Winn, J. N., Petigura, E. A., Morton, T. D., et al. 2017, AJ, 154, 270, doi: 10.3847/1538-3881/aa93e3
  * Yu & Tremaine (2001) Yu, Q., & Tremaine, S. 2001, AJ, 121, 1736, doi: 10.1086/319401
  * Zhao et al. (2022) Zhao, L. L., Kunovac, V., Brewer, J. M., et al. 2022, Nature Astronomy, doi: 10.1038/s41550-022-01837-2