# Aligning Planet-Hosting Binaries via Dissipative Precession in Circumstellar
Disks

Konstantin Gerbig Department of Astronomy, Yale University, New Haven, CT
06511, USA  Malena Rice Department of Astronomy, Yale University, New Haven,
CT 06511, USA  J.J. Zanazzi 51 Pegasi b Fellow  Astronomy Department and
Center for Integrative Planetary Science, University of California Berkeley,
Berkeley, CA 94720, USA  Sam Christian Department of Physics and Kavli
Institute for Astrophysics and Space Research, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  Andrew Vanderburg Department of Physics
and Kavli Institute for Astrophysics and Space Research, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA

###### Abstract

Recent observations have demonstrated that some subset of even moderately
wide-separation planet-hosting binaries are preferentially configured such
that planetary and binary orbits appear to lie within the same plane. In this
work, we explore dissipation during the protoplanetary disk phase, induced by
disk warping as the system is forced into nodal recession by an inclined
binary companion as a possible avenue of achieving orbit-orbit alignment. We
analytically model the coupled evolution of the disk angular momentum vector
and stellar spin vector under the influence of a distant binary companion. We
find that a population of systems with random initial orientations can appear
detectably more aligned after undergoing dissipative precession, and that this
process can simultaneously produce an obliquity distribution that is
consistent with observations. While dissipative precession proceeds
efficiently in close binaries, favorable system properties (e.g.,
$r_{out}\gtrsim 100$ AU, $\alpha\gtrsim 0.05$, and/or $M_{b}/M_{*}\gtrsim 1$)
are required to reproduce observed alignment trends at wider binary
separations $a_{\mathrm{b}}\gtrsim 450$ AU. Our framework further predicts
that circum-primary planets in systems with high stellar mass ratios should be
preferentially less aligned than planets in equal-mass stellar binary systems.
We discover tentative evidence for this trend in Gaia DR3 and TESS data. Our
findings suggest that dissipative precession may play a significant role in
sculpting orbital configurations in a sub-set of moderately-wide planet-
hosting binaries, but is likely not solely responsible for their observed
population-level alignment.

planetary alignment (1243), exoplanet dynamics (490), star-planet interactions
(2177), exoplanets (498), binary stars (154), exoplanet systems (484)

††software: NumPy (Harris et al., 2020), Matplotlib (Hunter, 2007), CMasher
(van der Velden, 2020)

##  1 Introduction

Many host stars of known exoplanets also possess a moderate- to wide-
separation stellar binary companion. The sky-projected orbital orientations
within such s-type binary systems are diverse, ranging from fully aligned
(Rice et al., 2023, 2024) to polar configurations with high stellar
obliquities (e.g. Winn et al., 2010; Sanchis-Ojeda et al., 2013; Esposito et
al., 2014; Rubenzahl et al., 2021; Albrecht et al., 2021; Hjorth et al.,
2021). The Gaia DR3 catalogue (Gaia Collaboration et al., 2022) allows stellar
orbit fitting of planet-hosting binaries to an unprecedented precision, and
has thus opened the pathway for population-level studies of their orbital
geometries. A key finding is the trend of planetary orbits to preferentially
appear aligned with the binary orbit (Christian et al., 2022; Dupuy et al.,
2022; Zhang et al., 2023; Lester et al., 2023; Rice et al., 2024; Christian et
al., 2024): a statistically significant preference that holds even up to a
binary separation of nearly $10^{3}$ AU. This trend can be interpreted as the
combined outcome of the systems’ birth orientations as resultant from star
formation (Bate et al., 2010; Bate, 2018) and secular interactions that may
occur ‘primordially’, between the binary companion and circumstellar disk
(Batygin, 2012; Batygin & Adams, 2013; Lai, 2014; Spalding et al., 2014;
Xiang-Gruess & Papaloizou, 2014; Zanazzi & Lai, 2017), and/or post planet
formation, between the binary companion and planet(s) (Naoz, 2016; Zhang et
al., 2018).

While star formation via disk fragmentation would naturally explain a
preferential emergence of aligned systems, it is unlikely to be the dominant
mode of binary formation outside of $\sim$ 500 AU (see e.g., Bonnell & Bate,
1994; Tobin et al., 2018; Kuruwita & Haugbølle, 2023; Offner et al., 2023). On
the other hand, core fragmentation of a (turbulent) molecular cloud core is
able to produce a much wider range of binary orbital separations (e.g., Padoan
& Nordlund, 2002), but results in largely randomly oriented initial
configurations (Bate et al., 2010).

We thus investigate the possibility of generating alignment — or rather,
damping misalignment — during the protoplanetary disk epoch. This is an
intuitive starting point, as gaseous protoplanetary disks possess viscosity
and thus offer a straight-forward avenue to dissipate orbital energy. Our work
is hence distinct from previous studies, such as Zhang et al. (2018), which
consider post-disk evolution of planetary embryos which may align (or anti-
align) with the accompanying binary orbit through chaotic effects – namely
close encounters or direct collisions. Such processes work within $\sim 200$
AU only (Zhang et al., 2018) and are thus likely insufficient to explain
observed alignment trends at wider separations.

The paradigm under scrutiny instead considers a protoplanetary disk that is
initially inclined with respect to the binary orbit. The binary companion
drives the disk into differential nodal recession (Batygin, 2012; Lai, 2014),
leading to warps. The resulting back-reaction torque evokes viscosity-
dependent damping terms (Foucart & Lai, 2014; Zanazzi & Lai, 2018a). At the
same time, the disk’s coupling to the stellar spin vector has been identified
as a potential source of obliquity excitation, which is expected if the
recession rate of the stellar spin about the disk falls below the disk’s
recession rate about the binary. Therefore, the same mechanism jointly sculpts
the spin-orbit distribution of exoplanets within binary star systems.

In this paper, we revisit and expand upon the considerations in Lai (2014) and
Zanazzi & Lai (2018a), where warp amplitudes are assumed to remain small such
that the disk is behaving almost like a rigid body. Unlike the individual
systems probed in these previous works, our analysis specifically seeks to
investigate the impact of dissipative precession mitigated by viscosity on a
population of disks. For this purpose, we co-evolve the stellar spin and disk
angular momentum vectors of a set of systems using a wide range of initial
orientations and track the expected distribution of sky-projected orbital
configuration angles.

In Section 2 we present the relevant equations of motion. Section 3 reviews
viscous dissipation in precessing disks and highlights our disk model within
the context of the current observational consensus. Section 4 introduces the
observer-oriented coordinate system, required to appropriately compare to
observations of sky-projected orientations. Section 5 contains the results of
our analysis. Section 6 considers observational data from Gaia DR3 and the
Transiting Exoplanet Survey Satellite (TESS), in the specific context of our
model’s predictions. Section 7 outlines additional effects which may sculpt
orbital configurations in planet-hosting binaries, and our conclusions can be
found in Section 8.

##  2 Nodal Recession due to Binary Potential

We consider a primary star with mass $M_{*}$, radius $R_{*}$, and rotation
rate $\Omega_{*}$ at the coordinate origin. The primary hosts a protoplanetary
disk and co-orbits with a stellar binary companion with mass $M_{\mathrm{b}}$
at a semi-major axis $a_{\mathrm{b}}$. The binary orbit, characterized by the
unit normal vector $\hat{\bm{l}}_{\mathrm{b}}$, torques the primary star’s
disk. This drives nodal recession of the disk and, if the primary star’s spin
axis is coupled to the disk, a corresponding recession of the stellar spin
vector.

To model this process, we follow Zanazzi & Lai (2018a) and track the evolution
of the unit normal vector of the disk orbital angular momentum
$\hat{\bm{l}}_{\mathrm{d}}$ and the unit vector of the primary star’s spin
axis $\hat{\bm{s}}$ via (see also Lai, 2014; Zanazzi & Lai, 2018b)

$\displaystyle\begin{split}\frac{\mathrm{d}\hat{\bm{l}}_{\mathrm{d}}}{\mathrm{d}t}=&-\tilde{\omega}_{\mathrm{ds}}(\hat{\bm{l}}_{\mathrm{d}}\cdot\hat{\bm{s}})\hat{\bm{s}}\times\hat{\bm{l}}_{\mathrm{d}}\\\
&-\tilde{\omega}_{\mathrm{db}}(\hat{\bm{l}}_{\mathrm{d}}\cdot\hat{\bm{l}}_{\mathrm{b}})\hat{\bm{l}}_{\mathrm{b}}\times\hat{\bm{l}}_{\mathrm{d}}+\left(\frac{\mathrm{d}\hat{\bm{l}}_{\mathrm{d}}}{\mathrm{d}t}\right)_{\mathrm{damp}},\end{split}$
(1) $\displaystyle\frac{\mathrm{d}\hat{\bm{s}}}{\mathrm{d}t}=$
$\displaystyle-\tilde{\omega}_{\mathrm{sd}}(\hat{\bm{s}}\cdot\hat{\bm{l}}_{\mathrm{d}})\hat{\bm{l}}_{\mathrm{d}}\times\hat{\bm{s}}+\left(\frac{\mathrm{d}\hat{\bm{s}}}{\mathrm{d}t}\right)_{\mathrm{damp}}.$
(2)

Here, $\tilde{\omega}_{\mathrm{ds}},\tilde{\omega}_{\mathrm{db}}$ and
$\tilde{\omega}_{\mathrm{sd}}$ are the recession rates of the disk around the
star, the disk around the binary, and the star around the disk, respectively,
and $({\mathrm{d}\hat{\bm{l}}_{\mathrm{d}}}/{\mathrm{d}t})_{\mathrm{damp}}$
and $({\mathrm{d}\hat{\bm{s}}}/{\mathrm{d}t})_{\mathrm{damp}}$ are damping
terms for the disk orbit normal and the primary spin, respectively. The binary
orbit is assumed be circular (eccentricity $e=0$), and
$\hat{\bm{l}}_{\mathrm{b}}$ is assumed to remain constant in time.

We define the mutual inclinations between the set of unit vectors as

$\displaystyle\theta_{\mathrm{sd}}$
$\displaystyle=\arccos(\hat{\bm{s}}\cdot\hat{\bm{l}}_{\mathrm{d}}),$ (3)
$\displaystyle\theta_{\mathrm{sb}}$
$\displaystyle=\arccos(\hat{\bm{s}}\cdot\hat{\bm{l}}_{\mathrm{b}}),$ (4)
$\displaystyle\theta_{\mathrm{db}}$
$\displaystyle=\arccos(\hat{\bm{l}}_{\mathrm{d}}\cdot\hat{\bm{l}}_{\mathrm{b}}).$
(5)

These angles track the relative star-disk ($\theta_{\mathrm{sd}}$), star-
binary ($\theta_{\mathrm{sb}}$), and disk-binary ($\theta_{\mathrm{db}}$)
evolution as the system evolves.

###  2.1 Recession rates

The binary-driven recession rates are dependent on the system’s stellar
separation, orientation, and disk properties. Following Zanazzi & Lai (2018b),
we consider a rigidly rotating protoplanetary disk with surface density
profile $\Sigma(r)\propto r^{-p}$, for $r$ between the inner and outer
truncation radii $r_{\mathrm{in}}$ and $r_{\mathrm{out}}$. The disk is assumed
to experience homologous mass loss, where the total disk mass evolves as (see
also Lai, 2014; Zanazzi & Lai, 2018b)

$\displaystyle M_{\mathrm{d}}=\frac{M_{\mathrm{d,0}}}{1+t/t_{\mathrm{v}}},$
(6)

with initial disk mass $M_{\mathrm{d,0}}$ and viscous timescale
$t_{\mathrm{v}}=r_{\mathrm{out}}^{2}/\nu$. It is useful to express kinematic
viscosity $\nu$ in terms of the Shakura-Sunyaev $\alpha$-parameter (Shakura &
Sunyaev, 1973) via $\nu=\alpha c_{\mathrm{s}}H$, where $c_{\mathrm{s}}$ is the
sound-speed, $H=c_{\mathrm{s}}/\Omega$ is the disk pressure scale height, and
$\Omega=\sqrt{GM_{*}/r^{3}}$ is the orbital frequency.

For this model, the precession frequencies are given in Zanazzi & Lai (2018b)
as

$\displaystyle\tilde{\omega}_{\mathrm{ds}}$
$\displaystyle=\frac{3(5/2-p)k_{\mathrm{q}}}{2(1+p)}\frac{R_{*}^{2}\bar{\Omega}_{*}^{2}}{r_{\mathrm{out}}^{1-p}r_{\mathrm{in}}^{1+p}}\sqrt{\frac{GM_{*}}{r_{\mathrm{out}}^{3}}},$
(7) $\displaystyle\tilde{\omega}_{\mathrm{db}}$
$\displaystyle=\frac{3(5/2-p)}{4(4-p)}\left(\frac{M_{\mathrm{b}}}{M_{*}}\right)\left(\frac{r_{\mathrm{out}}}{a_{\mathrm{b}}}\right)^{3}\sqrt{\frac{GM_{*}}{r^{3}_{\mathrm{out}}}},$
(8) $\displaystyle\tilde{\omega}_{\mathrm{sd}}$
$\displaystyle=\frac{3(2-p)k_{\mathrm{q}}}{2(1+p)k_{*}}\left(\frac{M_{\mathrm{d}}}{M_{*}}\right)\frac{\bar{\Omega}_{*}\sqrt{GM_{*}R_{*}^{3}}}{r_{\mathrm{out}}^{2-p}r_{\mathrm{in}}^{1+p}},$
(9)

where $k_{*}$ is the stellar spin normalization, $k_{\mathrm{q}}$ is the
stellar quadrupole moment normalization, and
$\bar{\Omega}_{*}=\Omega_{*}/\sqrt{GM_{*}/R_{*}^{3}}$ is the dimensionless
rotation rate of the primary. If the binary has nonzero eccentricity $e_{\rm
b}$, $a_{\rm b}$ in Eq. (8) is replaced with $a_{\rm b,eff}=a_{\rm
b}\sqrt{1-e_{\rm b}^{2}}$ (Anderson et al., 2016). Due to the mass loss of the
disk, $\tilde{\omega}_{\mathrm{sd}}$, initially the fastest precession rate,
falls below $\tilde{\omega}_{\mathrm{db}}$ within a few viscous timescales.
This implies a secular resonance is crossed ($\tilde{\omega}_{\rm
sd}\sim\tilde{\omega}_{\rm db}$) as the disk’s mass depletes, which excites
$\theta_{\mathrm{sd}}$ (and consequently stellar obliquities; see Lai 2014 for
discussion).

We can write the absolute value of the disk angular momentum as (Zanazzi &
Lai, 2018a)

$\displaystyle
L_{\mathrm{d}}=\frac{2-p}{5/2-p}M_{\mathrm{d}}\sqrt{GM_{*}r_{\mathrm{out}}}.$
(10)

The fact that the binary orbital angular momentum $L_{\mathrm{b}}$ satisfies
$L_{\mathrm{b}}\gg L_{\mathrm{d}}$ and also $L_{\mathrm{b}}\gg
S=k_{*}M_{*}R_{*}^{2}\Omega_{*}$ justifies treating the binary orbit as
static.

##  3 Viscous dissipation

In this work, we are primarily interested in understanding the joint behavior
of the spin-orbit angle $\theta_{\mathrm{sd}}$ and the orbit-orbit angle
$\theta_{\mathrm{db}}$ (as well as the corresponding sky-projected angles that
will be introduced in Section 4), the latter of which can only noticeably
evolve if the system has a means of dissipating energy. Indeed, any loss of
energy will move the system toward the lowest energy state – that is, full
alignment of the primary star’s spin axis, disk orbit, and binary orbit, as
seen in e.g. the Qatar-6 system (Rice et al., 2023, 2024). Intriguingly, a
polar configuration also has an absence of net-torque and thus constitutes an
unstable equilibrium.

In this section, we re-examine the ability of a precessing protoplanetary disk
to dissipate energy, quantifying the associated damping rates and dissipation
timescales.

###  3.1 Small amplitude disk warps

The formalism outlined in Sect. 2 assumes that the disk remains effectively
flat, in the sense that it precesses like a rigid body (also see Batygin &
Adams, 2013; Lai, 2014; Spalding & Batygin, 2015). This rigid body
approximation has to be enforced by effective communication between different
disk annuli through either self-gravity or hydrodynamical effects.

In this paper, we attribute this role to internal torques induced by small
warps in the disk profile – small amplitude bending waves (Papaloizou & Lin,
1995; Lubow & Ogilvie, 2000) – which dissipate energy when resisting the
external driving torque (Zanazzi & Lai, 2018a), and ultimately result in the
damping terms in Eqs. (1) and (2). In order for our model to remain self-
consistent, these warps have to remain small, which is the case as long as the
bending wave crossing time is shorter than disk precession timescales. Zanazzi
& Lai (2018a) show that this is indeed an appropriate limit for the adopted
framework. We thus proceed under the small warp approximation, which formally
presupposes that the time-dependent component of the warp unit vector, i.e.
the viscous twist, remains small: $|\bm{l}_{1}|\ll 1$. Note that this
eliminates the possibility for disk breaking, which otherwise can result in a
plethora of different dynamics (e.g., Larwood & Papaloizou, 1997; Lodato &
Price, 2010; Nixon & King, 2012; Doğan et al., 2015; Doǧan et al., 2018;
Smallwood et al., 2023).

Zanazzi & Lai (2018a) provide semi-analytic expressions for the damping terms
associated with small amplitude disk warps, i.e.

$\displaystyle\begin{split}\left(\frac{\mathrm{d}\hat{\bm{l}}_{\mathrm{d}}}{\mathrm{d}t}\right)_{\mathrm{visc}}&=\gamma_{\mathrm{b}}(\hat{\bm{l}}_{\mathrm{d}}\cdot\hat{\bm{l}}_{\mathrm{b}})^{3}\hat{\bm{l}}_{\mathrm{d}}\times(\hat{\bm{l}}_{\mathrm{b}}\times\hat{\bm{l}}_{\mathrm{d}})\\\
&+\gamma_{\mathrm{s}}(\hat{\bm{l}}_{\mathrm{d}}\cdot\hat{\bm{s}})^{3}\hat{\bm{l}}_{\mathrm{d}}\times(\hat{\bm{s}}\times\hat{\bm{l}}_{\mathrm{d}}),\end{split}$
(11)
$\displaystyle\left(\frac{\mathrm{d}\hat{\bm{s}}}{\mathrm{d}t}\right)_{\mathrm{visc}}$
$\displaystyle=-\frac{L_{\mathrm{d}}}{S}\gamma_{\mathrm{s}}(\hat{\bm{l}}_{\mathrm{d}}\cdot\hat{\bm{s}})^{2}\hat{\bm{s}}\times(\hat{\bm{s}}\times\hat{\bm{l}}_{\mathrm{d}}),$
(12)

where

$\displaystyle\begin{split}\gamma_{\mathrm{b}}=&1.26\cdot
10^{-9}\left(\frac{\alpha}{10^{-2}}\right)\left(\frac{h_{\mathrm{out}}}{0.1}\right)^{-2}\left(\frac{M_{\mathrm{b}}}{1M_{\odot}}\right)^{2}\\\
&\cdot\left(\frac{r_{\mathrm{out}}}{50\mathrm{AU}}\right)^{9/2}\left(\frac{a_{\mathrm{b}}}{300\mathrm{AU}}\right)^{-6}\left(\frac{M_{*}}{1M_{\odot}}\right)^{-3/2}\frac{2\pi}{\mathrm{yr}}\end{split}$
(13)

is the damping rate due to warps induced by the binary potential, and

$\displaystyle\begin{split}\gamma_{\mathrm{s}}=&2.04\cdot
10^{-10}\left(\frac{\alpha}{10^{-2}}\right)\left(\frac{h_{\mathrm{in}}}{0.1}\right)^{-2}\\\
&\cdot\left(1358\frac{r_{\mathrm{out}}}{r_{\mathrm{in}}}\right)^{p-1}\left(\frac{k_{\mathrm{q}}}{0.1}\right)^{2}\left(\frac{\bar{\Omega}_{*}}{0.1}\right)^{4}\left(\frac{M_{*}}{1M_{\odot}}\right)^{1/2}\\\
&\cdot\left(\frac{R_{*}}{2R_{\odot}}\right)^{4}\left(\frac{r_{\mathrm{in}}}{8R_{\odot}}\right)^{-4}\left(\frac{r_{\mathrm{out}}}{50\mathrm{AU}}\right)^{-3/2}\frac{2\pi}{\mathrm{yr}}\end{split}$
(14)

is the damping rate due to warps induced by an oblate star. Notably, we ignore
additional back-reaction terms, as the associated effect is dynamically
negligible compared to that of $\gamma_{\mathrm{b}}$ and $\gamma_{\mathrm{s}}$
(see Zanazzi & Lai, 2018a). In the above scaling, we write the disk aspect
ratio as $h=H/r$, such that at the inner and outer truncation radius,
$h_{\mathrm{in}}=H(r=r_{\mathrm{in}})/r_{\mathrm{in}}$ and
$h_{\mathrm{out}}=H(r=r_{\mathrm{out}})/r_{\mathrm{out}}$, respectively.

When plugging in fiducial model parameters (see the scaling relations in Eqs.
(13) and (14)), Zanazzi & Lai (2018a) find that, if $a_{\mathrm{b}}\gtrsim
200$ AU, viscous dissipation induced by disk warps is too weak to align the
disk with the binary orbit plane within typical disk lifetimes. However, their
modeled disks were relatively small and thus not necessarily representative of
a typical (early stage) protoplanetary disk in a binary, especially given the
high rate of large disks observed (e.g., Najita & Bergin, 2018; Boyden &
Eisner, 2020). In addition, larger effective viscosities than $\alpha=10^{-2}$
would further increase the damping rate. We discuss this prospect further in
the next section.

###  3.2 Disk damping with high viscosity

The viscosity of protoplanetary disks is relatively poorly constrained,
ranging from $\alpha\sim 10^{-4}$ from observations of edge-on protoplanetary
disks (Villenave et al., 2020; Rosotti, 2023) or numerical calculations (e.g.,
Flock et al., 2017) to $\alpha\sim 10^{-2}$, which well characterizes the
viscous evolution of average protoplanetary disks (Hartmann et al., 1998).
Indeed, Rafikov (2017) argues that this range is indicative of an intrinsic
viscosity distribution. Since these values were obtained by studying isolated
disks which are likely not warped, it is conceivable that the associated
$\alpha$-values are not necessarily appropriate for the disks in our study. In
particular, disk warps can have resonant interactions with inertial waves,
which result in a parametric instability (Gammie et al., 2000; Ogilvie &
Latter, 2013) that can increase viscosity and dissipation. If such an inertial
wave resonance successfully operates in the precessing and warped disks in our
study, the effective $\alpha$ values appearing in the damping rates in
Eqs.(13) and (14) may conceivably exceed the fiducial values in isolated disks
by orders of magnitude. Therefore, in the following analysis, we will consider
a viscosity of $\alpha=0.05$. For comparison, the fastest observed accretion
rates are also associated with theoretically puzzling high $\alpha$-values up
to 0.1 (Rafikov, 2017; Ansdell et al., 2018, see also Sect. 3.5).

The bending wave regime only holds if $\alpha<h$ (Papaloizou & Pringle, 1983).
If, on the other hand, $\alpha>h$, the disk lies in the diffusive regime. This
subsection shows that the disk damping rates in both regimes are identical, so
the same sets of evolutionary equations can be used regardless of the value of
$\alpha$ relative to $H/r$.

We start with the equations describing small warps in viscous disks, which
satisfy $|\kappa^{2}-\Omega^{2}|/\Omega^{2}\ll\alpha^{2}\ll 1$, where $\kappa$
is the epicyclic frequency. In this regime, the evolutionary equations for the
warp may be written as (Ogilvie, 1999; Dullemond et al., 2022)

$\displaystyle\Sigma r^{2}\Omega\frac{\partial{\bm{\hat{l}}}}{\partial t}$
$\displaystyle=\Sigma{\bm{T}}_{\rm
ext}+\frac{1}{r}\frac{\partial{\bm{G}}}{\partial r},$ (15)
$\displaystyle{\bm{G}}$ $\displaystyle\simeq\frac{\Sigma
H^{2}r^{3}\Omega^{2}}{4}\left(\frac{1}{\alpha}\frac{\partial{\bm{\hat{l}}}}{\partial
r}+\frac{3}{2}{\bm{\hat{l}}}{\bm{\times}}\frac{\partial{\bm{\hat{l}}}}{\partial
r}\right).$ (16)

Since $\alpha\ll 1$, we will approximate

${\bm{G}}\simeq{\bm{G}}_{\rm visc}=\frac{\Sigma
H^{2}r^{3}\Omega^{2}}{4\alpha}\frac{\partial{\bm{\hat{l}}}}{\partial r}.$ (17)

The pressure term ${\bm{G}}_{\rm press}=\frac{3\Sigma
H^{2}r^{3}\Omega^{2}}{8}{\bm{\hat{l}}}{\bm{\times}}\frac{\partial{\bm{\hat{l}}}}{\partial
r}$ adds slightly to the warp profile, but is otherwise uninteresting
dynamically. Taking

${\bm{\hat{l}}}(r,t)={\bm{\hat{l}}}_{\rm d}(t)+{\bm{l}}_{1}(r,t)$ (18)

with $|{\bm{l}}_{1}|\ll 1$, and adding external torques from the binary (Eq. 1
of Zanazzi & Lai, 2018a) and oblate star (Eq. 3 of Zanazzi & Lai, 2018a), the
same perturbation arguments can be repeated, substituting Equation (17) for
Eq. (17) of Zanazzi & Lai (2018a). This leads to a viscous twist of

${\bm{l}}_{1}=V_{\rm b}({\bm{\hat{l}}}_{\rm b}{\bm{\cdot}}{\bm{\hat{l}}}_{\rm
d}){\bm{\hat{l}}}_{\rm b}{\bm{\times}}{\bm{\hat{l}}}_{\rm d}+V_{\rm
s}({\bm{\hat{s}}}{\bm{\cdot}}{\bm{\hat{l}}}_{\rm
d}){\bm{\hat{s}}}{\bm{\times}}{\bm{\hat{l}}}_{\rm d}.$ (19)

Here, $V_{\rm b}$ is given by Eq. (31) of Zanazzi & Lai (2018a), while $V_{\rm
s}$ is given by Eq. (48) of Zanazzi & Lai (2018a). Comparing our Equation (19)
to Eq. (59) of Zanazzi & Lai (2018a), we see that the viscous twist in the
diffusive regime is identical to that in the bending wave regime. Hence,
Equations (11) and (12) can still be used when $\alpha\gtrsim H/r$, as long as
$|{\bm{l}}_{1}|\ll 1$.

Figure 1: Schematic of the studied configuration. The left panel shows a
protoplanetary disk with a binary companion. The three unit vectors
$\hat{\bm{s}},\hat{\bm{l}}_{\mathrm{d}}$ and $\hat{\bm{l}}_{\mathrm{b}}$
associated with stellar spin, disk orbital momentum, and binary orbital
momentum respectively — visible in all panels — are drawn in the observer
oriented coordinate system. The $z$-axis points towards the observer, and fix
the $y$-axis to lie parallel to $\hat{\bm{l}}_{\rm d}$, corresponding to
$i_{\rm d}=90^{\circ}$. Longitudes (shown in dashed lines in the right panel)
$\Omega_{\mathrm{s}}$ and $\Omega_{\mathrm{b}}$ define the angles between
$-\hat{\bm{z}}$ with the projections of $\hat{\bm{l}}_{\rm b}$ and
$\hat{\bm{s}}$ into the $\hat{\bm{x}}$-$\hat{\bm{z}}$ plane. Projected
inclinations (shown in dotted lines in the right panel) $\lambda_{\mathrm{s}}$
and $\lambda_{\mathrm{b}}$ define the angles between $\hat{\bm{l}}_{\rm d}$
with the projections of $\hat{\bm{l}}_{\rm b}$ and $\hat{\bm{s}}$ onto the
$\hat{\bm{x}}$-$\hat{\bm{y}}$ plane. Inclinations are displayed by solid arcs
in the center panel. The angles $\theta_{\rm db}$, $\theta_{\rm sd}$, and
$\theta_{\rm sb}$ define the angles between the three unit vectors
$\hat{\bm{l}}_{\rm d}$, $\hat{\bm{l}}_{\rm b}$, and $\hat{\bm{s}}$, while
$i_{\rm d}=90^{\circ}$, $i_{\rm b}$, and $i_{\rm s}$ define the angles between
the observer’s direction $\hat{\bm{z}}$ with $\hat{\bm{l}}_{\rm d}$,
$\hat{\bm{l}}_{\rm b}$, and $\hat{\bm{s}}$ respectively.

###  3.3 Suppression of Kozai-Lidov oscillations

If the binary inclination relative to the disk is sufficiently large – that
is, if it falls within the range $39^{\circ}<\theta_{\mathrm{db}}<141^{\circ}$
– the disk can undergo Kozai-Lidov oscillations (von Zeipel, 1910; Kozai,
1962; Lidov, 1962) where eccentricity and inclination of disk particles are
exchanged in a manner such that $\sqrt{1-e^{2}}\cos i$ remains constant at the
quadrupole level of approximation (Martin et al., 2014; Fu et al., 2015; Naoz,
2016; Lubow & Ogilvie, 2017). Using smoothed particle hydrodynamics
calculations, Smallwood et al. (2023) argue that highly inclined disks are
short-lived due to shocks that damp the disks’ relative tilt within a few
Kozai-Lidov timescales.

On the other hand, in protoplanetary disks, gas pressure forces slightly
eccentric orbits into apsidal precession, which suppresses Kozai-Lidov
oscillations (Wu & Murray, 2003). Zanazzi & Lai (2017) provide a criterion for
the suppression of Kozai-Lidov instability due to pressure-induced apsidal
precession, given by

$\displaystyle h_{\mathrm{out}}\gtrsim
0.08\cdot\left(\frac{5r_{\mathrm{out}}}{a_{\mathrm{b}}}\right)^{3/2}\left(\frac{M_{\mathrm{b}}}{M_{*}}\right)^{1/2}.$
(20)

Disk self-gravity also drives apsidal precession of disk annuli, likewise
working to suppress Kozai-Lidov oscillations (Batygin et al., 2011).
Specifically, if the disk mass exceeds (Zanazzi & Lai, 2017)

$\displaystyle M_{\mathrm{d}}\gtrsim
0.01M_{\mathrm{b}}\cdot\left(\frac{5r_{\mathrm{out}}}{a_{\mathrm{b}}}\right)^{3},$
(21)

Kozai-Lidov oscillations are suppressed.

In this work, we are primarily concerned with the evolution of massive disks
after the formation of the host star and the binary companion. Since we assume
homologous mass loss via Eq. (6), we can rephrase Eq. (21) into a criterion on
our time integration, i.e.

$\displaystyle t\lesssim
200t_{\mathrm{v}}\cdot\left(\frac{2M_{\mathrm{d,0}}}{M_{\mathrm{b}}}\right)\left(\frac{5r_{\mathrm{out}}}{a_{\mathrm{b}}}\right)^{3},$
(22)

which, for $t_{\mathrm{v}}=0.5$ Myr (see Sect. 3.5), lies well beyond the
expected lifetime of disks. As such, we solve the evolution of systems where
Eq. (21) holds such that Kozai-Lidov damping is suppressed by self-gravity. We
forgo including associated damping of relative inclinations in our
evolutionary equations. We will revisit the Kozai-Lidov mechanism in Section
7.2 when discussing post-planet-formation systems without a disk that acts as
suppressant.

###  3.4 Scaling relations

In order to reduce the number of free parameters in our model, we employ
scaling relations. Consider an equilibrium disk where the temperature scales
as (Chiang & Goldreich, 1997; Dullemond, 2000; Wu & Lithwick, 2021)

$\displaystyle T=140\mathrm{\
K}\cdot\left(\frac{M_{*}}{M_{\odot}}\right)^{-\frac{1}{7}}\left(\frac{L_{*}}{L_{\odot}}\right)^{\frac{2}{7}}\left(\frac{r}{1\mathrm{\
AU}}\right)^{-\frac{9}{20}}.$ (23)

Using an isothermal sound speed of $c_{\mathrm{s}}=\sqrt{k_{\mathrm{b}}T/\mu
m_{\mathrm{p}}}$, the aspect ratio at the disk outer edge is

$\displaystyle\begin{split}h_{\mathrm{out}}=0.08\cdot\left(\frac{M_{*}}{M_{\odot}}\right)^{-\frac{4}{7}}\left(\frac{L_{*}}{L_{\odot}}\right)^{\frac{1}{7}}\left(\frac{r_{\mathrm{out}}}{100\mathrm{AU}}\right)^{\frac{11}{40}},\end{split}$
(24)

assuming $\mu=2.33$ and provided the disk is locally isothermal, which is a
reasonable assumption for passively heated disks (Kratter & Murray-Clay,
2011). The disk aspect ratio at the inner edge using the same scaling is
$h_{\mathrm{in}}\sim 0.01$.

Hotter stars are also more massive. Specifically, we consider a mass-
luminosity relation $L_{*}\propto M_{*}^{A}$ with $A=3/2$ following Chachan &
Lee (2023), a choice motivated by MESA isochrone and stellar track models
(Dotter, 2016; Choi et al., 2016). We thus have

$\displaystyle
h_{\mathrm{out}}=0.08\cdot\left(\frac{M_{*}}{M_{\odot}}\right)^{-\frac{5}{14}}\left(\frac{r_{\mathrm{out}}}{100\mathrm{AU}}\right)^{\frac{11}{40}}.$
(25)

Note that, while more luminous stars host hotter disks, the corresponding
increase in stellar mass is associated with an increased vertical gravity that
leads to a net decrease in pressure scale height and consequently aspect ratio
at a given radius. Specifically, we can re-write Eq. (13) as

$\displaystyle\begin{split}\gamma_{\mathrm{b}}=&4.45\cdot
10^{-8}\left(\frac{\alpha}{10^{-2}}\right)\left(\frac{M_{\mathrm{b}}}{M_{*}}\right)^{2}\left(\frac{M_{*}}{1M_{\odot}}\right)^{\frac{17}{14}}\\\
&\cdot\left(\frac{r_{\mathrm{out}}}{100\mathrm{AU}}\right)^{\frac{79}{20}}\left(\frac{a_{\mathrm{b}}}{300\mathrm{AU}}\right)^{-6}\frac{2\pi}{\mathrm{yr}}.\end{split}$
(26)

For a constant binary mass ratio, having an overall massive and therefore hot
system is advantageous for strong damping, because the damping rate has a
stronger positive scaling with binary mass than its negative scaling with
primary mass. The novel insight of Eq. (26) is that inclination damping via
viscous dissipation operates more efficiently when the companion star is more
massive than the host star. In other words, if the mechanism studied in this
paper is the reason for observed alignment trends in planet-hosting binaries
found in Christian et al. (2022); Dupuy et al. (2022); Rice et al. (2024),
then circum-secondary planets (i.e. those orbiting the less massive binary
component) are expected to have experienced more efficient alignment than
circum-primary planets. We will test this hypothesis in Sect. 5.2.

###  3.5 Disk model in the context of observed protoplanetary disks

Our scaling for $h_{\rm out}$ in Eq. (25) also sets the disk viscous time
$t_{\rm v}\sim r^{2}/\nu$ evaluated at $r_{\mathrm{out}}$, i.e.,

$\displaystyle t_{\rm v}\sim 0.50\ {\rm
Myr}\cdot\left(\frac{0.05}{\alpha}\right)\left(\frac{M_{*}}{M_{\odot}}\right)^{\frac{3}{14}}\left(\frac{r_{\rm
out}}{100\ {\rm AU}}\right)^{\frac{19}{20}}.$ (27)

Because we take $r_{\rm out}=100\ {\rm AU}$ for all of our disks, and $t_{\rm
v}$ scales weakly with $M_{\star}$, we approximate $t_{\rm v}\approx{\rm
constant}=0.5\ {\rm Myr}$ for the rest of this work.

The initial mass loss of a protoplanetary disk is primarily set by its viscous
evolution, i.e. Eq. (6). On the other hand, its ultimate dispersal tends to be
set by the photoevaporative mass loss rate, where far and extreme ultraviolet
photons from the disk’s host star produce hot electrons via photochemical
reactions (e.g., Ercolano et al., 2009; Owen et al., 2011). Subsequent
thermalization heats the gas to a degree sufficient to escape the star’s
potential well. Since X-ray luminosities increase with stellar mass, so does
the photoevaporation timescale. Specifically, Komaki et al. (2021) find a mass
loss rate of $\dot{M}\propto M_{*}^{2}$. We neglect this scaling in this paper
and instead assume a mass-independent disk lifetime of $t_{\mathrm{l}}=10^{7}$
yr. If considered, this dependence would influence the expectations for
alignment distribution such that more massive systems would damp more quickly
but for a shorter period of time, whereas less massive systems would damp
inclinations more slowly but for a longer period of time. Combined with the
scaling relationship $\gamma_{\mathrm{b}}\propto M_{*}^{17/14}$ between
damping rate and total stellar mass in Eq. (26), we expect this effect to be
of negligible importance. We thus focus our study on the binary mass ratio —
on which the damping rate depends much more strongly — rather than the total
stellar mass.

Table 1: Overview of constant parameters in our model. Other prescribed parameters are the masses of central star $M_{*}$ and binary companion $M_{\mathrm{b}}$, as well as their separation $a_{\mathrm{b}}$. The aspect ratio at the inner and outer truncation radii, $h_{\mathrm{in}}$ and $h_{\mathrm{out}}$, are calculated via Eq. (25). Symbol | Meaning | Value
---|---|---
$R_{\mathrm{s}}$ | central star radius | $2R_{\odot}$
$\overline{\Omega}_{\mathrm{*}}$ | normalized central star rotation rate | $0.1$
$M_{\mathrm{d,0}}$ | initial disk mass | $0.3M_{*}$
$r_{\mathrm{in}}$ | inner truncation radius |  $0.1$ AU
$r_{\mathrm{out}}$ | outer truncation radius |  $100$ AU
$p$ |  $\Sigma$ power law index | $1.0$
$k_{\mathrm{q}}$ | quadrupole coefficient | 0.1
$k_{*}$ | spin normalization coefficient | 0.2
$t_{\mathrm{v}}$ | viscous timescale |  $0.5$ Myr
$\alpha$ | disk viscosity | $0.05$
$t_{\mathrm{l}}$ | disk lifetime |  $10$ Myr
Figure 2: Mass accretion rate and total disk mass for observed Class II disks
in Manara et al. (2023), with triangles corresponding to upper bounds in disk
mass. Total disk mass is calculated multiplying the recorded mm continuum mass
by a gas-to-dust mass ratio of 100. In color, we plot the trajectory of our
disk model in this space following the viscous mass loss prescription in Eq.
(6) and the fiducial parameters in Tab. 1. For the model, accretion rate
$M_{\mathrm{acc}}=|\mathrm{d}M_{\mathrm{d}}/\mathrm{d}t|$ is assumed to be due
to viscous evolution only, as opposed to wind-driven accretion. This is not
necessarily the case for the observed disks (see Manara et al., 2023).

Table 1 shows an overview of constant parameters in our model, largely adopted
from Zanazzi & Lai (2018a). Notable changes are the increased disk viscosity
(see Sect. 3.2) and a relatively large outer truncation radius of
$r_{\mathrm{out}}=100$ AU. This value remains consistent with observations,
which show a large variety of disk sizes from tens to hundreds of AU (Najita &
Bergin, 2018). The largest dust disks may extend up to $\sim 200$ AU (Hendler
et al., 2020), and the gas disk is thought to be a factor of $\sim 2-3$ larger
(Ansdell et al., 2018). For such large disks, a viscosity parameter of
$\alpha\sim 3\cdot 10^{-3}$ would, per Eq. (26), yield the same damping rate
as for our fiducial model with $\alpha=0.05$. The choice of
$r_{\mathrm{out}}=100\mathrm{AU}$ also ensures that all models considered in
this paper strictly adhere to $a_{\mathrm{b}}/r_{\mathrm{out}}\geq 4$,
avoiding disk truncation by the binary orbit. Note that, if the binary’s orbit
were eccentric, $a_{\rm b,eff}/r_{\rm out}=a_{\rm b}\sqrt{1-e_{\rm
b}^{2}}/r_{\rm out}\geq 4$ would avoid tidal truncation until $e_{\rm
b}\gtrsim 0.7$, since the disk’s tidal truncation radius $r_{\rm tide}\approx
0.35(1-e_{\rm b})a_{\rm b}=0.35[(1-e_{\rm b})/(1+e_{\rm b})]^{1/2}a_{\rm
b,eff}$ (e.g. Miranda & Lai, 2015).

The initial stages of the calculated evolution aim to reflect the transition
from Class I to Class II, hence the high initial disk mass of
$M_{\mathrm{d,0}}=0.3M_{*}$ (compare with e.g. Andersen et al., 2019), which
is then rapidly reduced following the utilized prescription for homologous
mass loss. This disk mass renders the outer regions of the modeled disk
(beyond $45$ AU) gravitationally unstable for $\lesssim 0.8t_{\mathrm{v}}\sim
0.4$ Myr, possibly increasing turbulence $\alpha$ (e.g., Gammie, 2001; Booth &
Clarke, 2019) and thus the damping rate $\gamma_{\mathrm{b}}$. We show the
trajectory of our model in total disk mass $M_{\mathrm{d}}$ – accretion rate
$M_{\mathrm{acc}}$ space in Fig. 2. For comparison, we show observational data
of Class II disks with corresponding measurements compiled by Manara et al.
(2023) (their Fig. 7). Disk (gas) mass $M_{\mathrm{d}}$ is calculated from the
dust mass traced by mm continuum emission assuming an interstellar medium-like
solid-to-gas ratio of $1\%$.

This approach relies on assumptions on uncertain dust opacities (e.g., Wright,
1987; Andrews & Williams, 2005; Birnstiel et al., 2018; Liu et al., 2022;
Miotello et al., 2023), and it presupposes the mass distribution of solids to
peak at mm sizes (Williams & Cieza, 2011; Manara et al., 2023). If, however,
radial drift, or aggregation via dust growth, and/or planet formation operate
even moderately fast, some fraction of the total solid mass may be either lost
or hidden in larger dust grains or planetesimals (Hughes & Armitage, 2012;
Gerbig et al., 2019), leading to very different true solid-to-gas mass ratio
and a potentially substantial under-estimation of total disk masses (e.g.,
Birnstiel et al., 2012; Powell et al., 2019). With this in mind, Fig. 2
characterizes our model as a moderately massive and large, relatively slowly
accreting disk that is well within the bounds of plausability established by
the current observational and theoretical consensus.

On first glance, the low accretion rates of the model relative to typical
disks in the corresponding mass range appear to be inconsistent with the
fairly high viscosity of $\alpha=0.05$. Indeed, naively taking for example Eq.
(4) in Rafikov (2017), i.e.

$\displaystyle\alpha\sim\frac{M_{\mathrm{acc}}}{M_{\mathrm{d}}}\frac{\Omega
r^{2}}{c_{\mathrm{s}}^{2}}$ (28)

evaluated at $r_{\mathrm{out}}$, one would infer an $\alpha$ that decreases in
time from $\sim 0.08$ to $\sim 0.005$. This is because our model assumes a
polynomial mass loss via Eq. (6) as opposed to the exponential decay
assumption that underlies Eq. (28) – a reminder that subtleties in the assumed
viscous evolution model may lead to appreciable differences when connecting
accretion rates and viscosity.

To conclude this section, in the construction of our fiducial model, we have
deliberately selected parameters that lead to favorable conditions for the
studied mechanism. The parameter choices are plausible for a sub-set of disks,
but should not be taken as representative for the entire population. This will
allow us to scrutinize the impact of viscous dissipation during binary-driven
disk precession on inclination and obliquity distributions. In Sect. 8, we
will comment on the implications of this choice on the mechanism’s ability to
explain observed alignment trends in exoplanet-hosting binaries.

##  4 Observer-oriented coordinate system

Figure 3: Possible trajectories in $i_{\mathrm{b}}$-$\lambda_{\mathrm{s}}$
space for sets of $\theta_{\mathrm{db}}$, $\theta_{\mathrm{sd}}$ as the
observer’s longitude relative to the stellar spin $\Omega_{\mathrm{s}}$
changes (see Fig. 1). From left to right, the four panels depict the following
configurations: close-to-full alignment, close-to-polar binary and moderately
misaligned spin, moderately misaligned binary and close-to-polar spin, and,
close-to-polar binary and retrograde close-to-polar spin. Each trajectory is
mapped out by 100 points spaced evenly from $\Omega_{\mathrm{s}}=0$ to
$\Omega_{\mathrm{s}}=\pi$. Their density thus can be interpreted as the
probability $p(i_{\mathrm{b}},\lambda_{\mathrm{s}})$. The true angles
$\theta_{\mathrm{db}}$, $\theta_{\mathrm{sd}}$ define a forbidden region (gray
shaded) which can never been observed. For given $\theta_{\mathrm{db}}$,
$\theta_{\mathrm{sd}}$, the third invariant angle $\theta_{\mathrm{sb}}$ is
geometrically constrained via
$\mathrm{min}(\theta_{\mathrm{sb}})=|\theta_{\mathrm{db}}-\theta_{\mathrm{sd}}|\leq\theta_{\mathrm{sb}}\leq\pi-|\pi-(\theta_{\mathrm{db}}+\theta_{\mathrm{sd}})|=\mathrm{max}(\theta_{\mathrm{sb}})$.
Each panel shows two trajectories corresponding to
$\mathrm{max}(\theta_{\mathrm{sb}})$ (circles) and
$\langle\theta_{\mathrm{sb}}\rangle=[\mathrm{max}(\theta_{\mathrm{sb}})-\mathrm{min}(\theta_{\mathrm{sb}})]/2$
(crosses). Note that the latter is an average of geometrically possible
$\theta_{\mathrm{sb}}$ defined for visualization only and is unrelated to
physical effects.

To directly compare the results of this work with observations, it is useful
to introduce an ‘observer-oriented’ coordinate system (compare to Fabrycky &
Winn, 2009). In this frame of reference, the coordinate origin is at the
location of the primary star, and the $\hat{\bm{z}}$ vector points towards the
observer, so that $\hat{\bm{x}}-\hat{\bm{y}}$ define the sky plane. We define
the line-of-sight inclinations of the three unit vectors as

$\displaystyle i_{\mathrm{d}}$
$\displaystyle=\arccos(\hat{\bm{l}}_{\mathrm{d}}\cdot\hat{\bm{z}}),$ (29)
$\displaystyle i_{\mathrm{s}}$
$\displaystyle=\arccos(\hat{\bm{s}}\cdot\hat{\bm{z}}),$ (30) $\displaystyle
i_{\mathrm{b}}$
$\displaystyle=\arccos(\hat{\bm{l}}_{\mathrm{b}}\cdot\hat{\bm{z}}).$ (31)

Figure 1 shows the set of unit vectors in this coordinate system. Because the
transiting planet forms in the disk, whose inclination we assume to evolve
negligibly after disk dispersal, we fix $i_{\rm d}=90^{\circ}$, and are free
to set the orientation of the sky plane so that the $y$-axis is parallel to
$\hat{\bm{l}}_{\rm d}$. This defines the observer-oriented reference plane we
will use from here on. Thus, $\theta_{\mathrm{sd}}$ is the planet’s true
obliquity (conventionally denoted as $\psi$), and $\theta_{\mathrm{db}}$
characterizes the ‘orbit-orbit’ angle, i.e. the true inclination between
planetary orbit and binary orbit.

Fabrycky & Winn (2009) provide the angle transformations, i.e.

$\displaystyle\tan\lambda_{\mathrm{s}}$
$\displaystyle=\frac{\sin\theta_{\mathrm{sd}}\sin\Omega_{\mathrm{s}}}{\cos\theta_{\mathrm{sd}}\sin
i_{\mathrm{d}}+\sin\theta_{\mathrm{sd}}\cos\Omega_{\mathrm{s}}\cos
i_{\mathrm{d}}},$ (32) $\displaystyle\tan\lambda_{\mathrm{b}}$
$\displaystyle=\frac{\sin\theta_{\mathrm{db}}\sin\Omega_{\mathrm{b}}}{\cos\theta_{\mathrm{db}}\sin
i_{\mathrm{d}}+\sin\theta_{\mathrm{db}}\cos\Omega_{\mathrm{b}}\cos
i_{\mathrm{d}}},$ (33)

with $0\leq\lambda_{\mathrm{s}},\lambda_{\mathrm{b}}\leq\pi$. Here,
$\Omega_{\mathrm{s}}$ and $\Omega_{\mathrm{b}}$ are the longitudinal angles of
stellar spin vector and binary orbit normal (see Fig. 1). For
$i_{\mathrm{d}}=90^{\circ}$,
$\tan\lambda_{\mathrm{s}}=\tan\theta_{\mathrm{sd}}\Omega_{\mathrm{s}}$ and
$\tan\lambda_{\mathrm{b}}=\tan\theta_{\mathrm{db}}\Omega_{\mathrm{b}}$

We also introduce an additional longitudinal angle, defined as the angle
between $\hat{\bm{s}}$ and $\hat{\bm{l}}_{\rm b}$ after projection onto the
$\hat{\bm{x}}$-$\hat{\bm{z}}$ plane:

$\displaystyle\Omega_{\mathrm{b}}=\Omega_{\mathrm{s}}+\Delta\Omega,$ (34)

It can be shown that (see Appendix A),

$\displaystyle\cos\Delta\Omega=\frac{\cos\theta_{\mathrm{sb}}}{\sin\theta_{\mathrm{db}}\sin\theta_{\mathrm{sd}}}-\frac{1}{\tan\theta_{\mathrm{db}}\tan\theta_{\mathrm{sd}}}.$
(35)

We will take the azimuthal orientation of our system to be random uniform,
which corresponds to $p(\Omega_{\mathrm{s}})=\pi^{-1}$. As the observables of
interest are $\lambda_{\mathrm{s}}$ and $i_{\mathrm{b}}$ (see Rice et al.,
2024), we are interested in computing the probability distributions
$p(\lambda_{\mathrm{s}}|\theta_{\mathrm{sd}},\Omega_{\mathrm{s}})$, and
$p(i_{\mathrm{b}}|\theta_{\mathrm{sb}},\theta_{\mathrm{sd}},\theta_{\mathrm{db}},\Omega_{\mathrm{s}})$
or $p(\cos
i_{\mathrm{b}}|\theta_{\mathrm{sb}},\theta_{\mathrm{sd}},\theta_{\mathrm{db}},\Omega_{\mathrm{s}})$.
The binary orbit inclination $i_{\mathrm{b}}$ relates to
$\lambda_{\mathrm{b}}$ via (Fabrycky & Winn, 2009)

$\displaystyle\sin
i_{\mathrm{b}}=\frac{\sin\theta_{\mathrm{db}}\sin\Omega_{\mathrm{b}}}{\sin\lambda_{\mathrm{b}}},$
(36)

while the inclination of the spin vector relative to the sky plane is

$\displaystyle\sin
i_{\mathrm{s}}=\frac{\sin\theta_{\mathrm{sd}}\sin\Omega_{\mathrm{s}}}{\sin\lambda_{\mathrm{s}}}.$
(37)

Figure 3 visualizes possible trajectories in
$i_{\mathrm{b}}$-$\lambda_{\mathrm{s}}$ space (observed binary inclination and
stellar obliquity), for four sets of $\theta_{\mathrm{db}}$ and
$\theta_{\mathrm{sd}}$ (‘true’ binary inclination and stellar obliquity). For
all $\Omega_{\mathrm{s}}$, $i_{\mathrm{b}}\geq\pi/2-\theta_{\mathrm{db}}$. For
stellar spins in prograde rotation relative to $\hat{\bm{l}}_{\mathrm{d}}$,
i.e. if $\theta_{\mathrm{sd}}<\pi/2$, then
$\lambda_{\mathrm{s}}\leq\theta_{\mathrm{sd}}$ always. For retrograde stellar
spins with $\theta_{\mathrm{sd}}>\pi/2$,
$\lambda_{\mathrm{s}}\geq\theta_{\mathrm{sd}}$.

###  4.1 Initial configuration

The three main mechanisms for forming stellar binary systems are disk
fragmentation, turbulent fragmentation, and dynamical capture, which are
associated with different expectations for the primordial inclination between
binary orbit and stellar spin orbit.While dynamical capture results in largely
random inclinations, disk fragmentation is expected to initialize the true
binary inclination as preferentially aligned. Since those two mechanisms are
typically invoked for the formation of very wide ($>10^{4}$ AU) (Tokovinin,
2017) and relatively close ($<500$ AU) (Krumholz et al., 2007; Tobin et al.,
2018) binaries, respectively, we instead focus on turbulent fragmentation,
which is shown to manufacture binaries with orbital separations throughout our
range of interest, i.e. $400<a_{\mathrm{b}}<1000$ AU (Padoan & Nordlund, 2002;
Bate, 2018; Lee et al., 2017), and with, to zeroth order, random orientations.

We thus assume the orientation of stellar spin angular momentum vector to be
random uniform. The corresponding probability density function for the angle
is given by $p[\theta_{\mathrm{db}}(t=0)]\propto\sin{|\theta_{\mathrm{db}}|}$.
Note that, in our model, retrograde binary orbits behave identically to
prograde orbits. We will focus our calculations on initial misalignments
$0^{\circ}\leq\theta_{\mathrm{db}}\leq 90^{\circ}$.

###  4.2 Example model

Figure 4: Evolution of disk-binary angle $\theta_{\mathrm{db}}$ (left panel)
and stellar spin- disk angle $\theta_{\mathrm{sd}}$ (center panel) for models
initialized at different inclinations. The probability map
$p(i_{\mathrm{b}},\lambda_{\mathrm{s}})$ is shown in the right panel (see
Sect. 4). The depicted models use our fiducial parameter set for an equal-mass
binary $M_{*}=M_{\mathrm{b}}=M_{\odot}$ with a binary separation of
$a_{\mathrm{b}}=600$ AU (see Tab.1).

We calculate the evolution of the vectors $\hat{\bm{l}}_{\mathrm{d}}$ and
$\hat{\bm{s}}$ over $t_{\mathrm{l}}=10^{7}$ yr, a fiducial value for the
lifetime of the gas disk consistent with expectations from photoevaporation
models (Ercolano et al., 2009; Owen et al., 2012). We then take the final
values of $\theta_{\rm sd}$, $\theta_{\rm sb}$, and $\theta_{\rm db}$ as
inputs for Eqs. (32) and (36) to calculate probability density functions for
the angles $\lambda_{\mathrm{s}}$ and $i_{\mathrm{b}}$. Note, that resonance
crossing can put the system into a state of oscillating $\theta_{\mathrm{sd}}$
with amplitudes up to tens of degrees (e.g., Lai, 2014). The phase of this
oscillation at $t=t_{\mathrm{l}}$ and the associated final value for the spin
orbit angle are sensitive to the initial condition, requiring sufficiently
fine sampling of initial binary inclinations $\theta_{\mathrm{db}}$.

Fig. 4 shows the evolution of the invariant angles $\theta_{\mathrm{db}}$ and
$\theta_{\mathrm{sd}}$ (left and center panel respectively), where color
corresponds to the initial $\theta_{\mathrm{db}}$, for our set of default
parameters depicted in Table 1, as well as $M_{*}=M_{\mathrm{b}}=M_{\odot}$,
at a separation of $a_{\mathrm{b}}=600$ AU. As studied in Zanazzi & Lai
(2018a, b), the depletion of disk material causes
$\tilde{\omega}_{\mathrm{sd}}$ to fall below $\tilde{\omega}_{\mathrm{db}}$,
triggering a resonance crossing. The right panel of Fig. 4 shows the
probability density map $p(i_{\mathrm{b}},\lambda_{\mathrm{s}})$ in arbitrary
units. Most initial alignments lead obliquity excitation into a state of
oscillating $\theta_{\mathrm{sd}}$. Many systems end up with a retrograde
spin-orbit angle $\theta_{\mathrm{sd}}>\pi/2$, which results in observed
trajectories in $i_{\mathrm{b}}$-$\lambda_{\mathrm{s}}$\- space with
$\pi>\lambda_{\mathrm{s}}>\pi/2$ (see Fig. 3). Systems not experiencing
significant excitation of obliquity are mostly those initialized with close-to
polar $\theta_{\mathrm{db}}$. These correspond to
$i_{\mathrm{b}}$-$\lambda_{\mathrm{s}}$\- space trajectories similar to the
one shown in the second from the left panel in Fig. 3, thus populating the
bottom part of the probability density map
$p(i_{\mathrm{b}},\lambda_{\mathrm{s}})$ in the right panel of Fig. 4. While
we recover the result in Zanazzi & Lai (2018a) – namely, that viscous damping
for the chosen set of parameters is too weak to fully drive individual disks
towards alignment with the binary plane – the population of disks as a whole
shifts away from an isotropic distribution towards a more aligned
distribution. We will elaborate on this finding further in the following
section.

##  5 Results

Figure 5: Probability of measuring $\cos(i_{\mathrm{b}})$ for binary
separations between 400 AU and 800 AU. Different panels correspond to
different binary configurations: circum-primary disk with less massive stellar
companion (left), equal mass binary (center), and circum-secondary disk with
more massive stellar companion (right). Figure 6: Same as Fig. 5 but for sky-
projected obliquity $\cos(\lambda_{\mathrm{s}})$. The predicted distribution,
specifically for equal-mass-ratio binaries, is broadly consistent with
observed trends in Dong & Foreman-Mackey (2023); Siegel et al. (2023) and Rice
et al. (2024). Figure 7: Cumulative probability map of equal-mass-ratio models
shown in the center panels of Figs. 5 and 6. In contrast to the right panel of
Fig. 4, this probability map covers $0^{\circ}\leq i_{\mathrm{b}}\leq
180^{\circ}$, achieved by mirroring at $i_{\mathrm{b}}=90^{\circ}$. We also
set the model prediction side-by-side with the data of the ‘strict’ sample of
planet-hosting binaries with obliquity measurements in Rice et al. (2024).
White dots belong to systems with binary separations that fall within
$400\,\mathrm{AU}<a_{\mathrm{b}}<800\,\mathrm{AU}$, whereas green dots are
systems either inside or outside this range.

We aim to characterize the maximum binary separation at which dissipative
precession significantly causes deviations from isotropy in the final
distribution of line-of-sight binary inclination $i_{\mathrm{b}})$ and sky-
projected stellar obliquity $\lambda_{\mathrm{s}}$, as well as the expected
variations for circum-primary versus circum-secondary planets in binary
systems with different mass ratios. For this purpose, Figs. 5 and 6 show the
probabilities of measuring $\cos(i_{\mathrm{b}})$ and $\lambda_{\mathrm{s}}$,
respectively, for three binary mass configurations across binary separations
ranging from 400-800 AU.

###  5.1 Equal-mass binary

The damping rate dramatically weakens with binary separation, such that for
typical disk parameters, even moderately wide binaries have damping timescales
that exceed expected disk lifetimes (see Fig. 5). Nonetheless, a population of
disks subject to this ‘weak’ damping can still have a distribution of sky-
projected angles that deviate from random orientations. To quantify the
maximum binary separation for which we reasonably expect this to occur, we
first consider equal-mass binaries with $M_{*}=M_{\mathrm{b}}=M_{\odot}$.

The center panel of Fig. 5 demonstrates that, for the fiducial setup, up until
a binary separation of $\sim$ 500 AU our model predicts an orientation
distribution that deviates from isotropy – specifically, an excess of aligned
systems and a lack of misaligned systems. Despite parameter choices that are
relatively favorable for alignment, this cut-off is short of the boundary
between preferentially aligned and randomly oriented systems, which Christian
et al. (2022) reported to be located between 600 and 800 AU. Moreover, in our
study, the inclinations of systems initialized in a near-polar configuration
($\cos i_{\rm b}\sim 1.0$) remain largely unaffected by the inclination
damping studied in this paper. While this feature has not been identified in
observations (Christian et al., 2022; Dupuy et al., 2022), where the dearth of
misaligned systems may persist at polar orientations, it may be a useful
indicator for efficient dissipation if seen in future observations that
incorporate larger sample sizes.

For these reasons, the herein studied mechanism alone can likely not explain
the observed orientations of planet-hosting binaries, but instead needs to be
accompanied by an additional process such as primordial alignment of the disk
or post planet formation reorganization. We expand on this in Sect. 7.

Obliquity excitation is less sensitive to semi-major axis (center panel of
Fig. 6) and still occurs up to $a_{\mathrm{b}}=800$ AU, although for a smaller
fraction of systems than at smaller binary separations. At larger binary
separations, secular resonance ($\tilde{\omega}_{\rm
sd}\sim\tilde{\omega}_{\rm db}$) occurs at later times, prolonging the
excitation of $\lambda_{\rm s}$ to timescales longer than those given by Eq.
(22). For sufficiently large $a_{\rm b}$, the binary’s torque becomes too weak
to excite $\lambda_{\rm s}$. This is why the sky projected obliquity
$\lambda_{\rm s}$ tends to decrease – or $\cos\lambda_{\rm s}$ tends to
increase – as $a_{\rm b}$ gets larger.

Overall, the predicted distribution qualitatively agrees with recently
reported obliquity distributions, both specifically for planet-hosting
binaries (Rice et al., 2024) and also across the full census of systems with
spin-orbit measurements (Albrecht et al., 2021; Dong & Foreman-Mackey, 2023;
Siegel et al., 2023). For large binary separations, most systems do not
experience obliquity excitation, and the stellar spin and disk orbit normal
remain aligned. Misaligned systems are expected for a wider range of initial
orientations at smaller binary separations. Albrecht et al. (2021) report an
overabundance of systems that are perpendicular ($\cos\lambda_{s}\sim 0$) in
their spin-orbit orientations (although see also Siegel et al. (2023) and Dong
& Foreman-Mackey (2023), which show that this peak is not significant in a
hierarchical Bayesian framework). A broad peak near $\cos\lambda_{s}=0$ may be
consistent with our model if resonant crossing preferentially excites
obliquities to near-polar configurations (see also Vick et al., 2023). Indeed,
the systems with significant spin-orbit misalignment tend to be in a state of
high-amplitude oscillation of $\theta_{\mathrm{sd}}$ around a near-polar value
(see Fig. 4). The herein utilized parameter set, however, does not damp these
oscillations to the mid-point value sufficiently fast, ultimately resulting in
a more spread-out distribution of $\cos(\lambda_{\mathrm{s}})$ in Fig. 6.

The cumulative probabilities $p(i_{\mathrm{b}},\lambda_{\mathrm{s}})$ of all
models with equal stellar mass ratio are shown in Fig. 7 in comparison to the
observed joint obliquity – binary inclination – distribution reported in Rice
et al. (2024). For this purpose, Fig. 7 mirrors the probability map at
$i_{\mathrm{b}}=90^{\circ}$. We only show the 40 systems in the ‘strict’
sample from Rice et al. (2024), out of which 8 have binary separation within
400 and 800 AU (red dots). The fiducial model and the (sparse) data have the
following trends in common: most systems appear close to spin-orbit alignment,
including a pile-up of fully aligned systems with a wide range of allowed
obliquities. We are not concerned that the shape of the probability map does
not accurately track the data points from Rice et al. (2024), as (1) the data
of planet-hosting binaries with obliquity measurements are sparse, (2)
obliquity $\lambda_{\mathrm{s}}$ (arguably unlike binary inclination
$i_{\mathrm{b}}$) can be affected by a wide variety of other dynamical effects
that operate during or after the protoplanetary disk phase (e.g., Anderson &
Lai, 2018; Millholland & Batygin, 2019; Li, 2021; Vick et al., 2023), and (3)
the model was (intentionally) not set up as a population synthesis, but to
highlight potential observable implications of the physics of dissipative
precession in the circumstellar disk. Such a population synthesis would need
to carefully consider the types of systems that underlie the observations,
explore a larger parameter space, and involve some treatment of post-disk
dynamical mechanisms that affect observables. We defer this to future work.

###  5.2 Hierarchical binary

We showed that disks in equal-mass-ratio binaries, on a population level, can
noticeably experience alignment. Yet, tension remains: the viscous dissipation
mechanism, even with optimistic model parameters, becomes inefficient for
exceedingly wide binaries ($a_{\mathrm{b}}\gtrsim 600$ AU). Consequently, the
study so far is an insufficient argument for alignment of exoplanet-hosting
binaries via viscous dissipation during the binary-driven recession of
protoplanetary disks, and other potential observable consequences must be
examined. The mass ratio dependence of Eq.(26) provides this opportunity.

We consider a hierarchical binary: specifically, a system with mass ratio of
$M_{1}/M_{2}=4$. Figures 5 and 6 show corresponding orientation distributions
in inclination and obliquity for the circum-primary systems (disk host star
mass $M_{*}=2M_{\odot}$ and companion star mass $M_{\mathrm{b}}=0.5M_{\odot}$)
in the left panels, and for circum-secondary systems (disk host star mass
$M_{*}=0.5M_{\odot}$ and companion star mass $M_{\mathrm{b}}=2M_{\odot}$) in
the right panels.

Our key finding from this exercise is that, for circum-primary systems, the
final inclination distribution more closely resembles an isotropic
distribution than for equal-mass-ratio binaries, even at small projected
separations of $a_{\mathrm{b}}\sim 400$ AU. For circum-secondary systems, by
contrast, the inclination distribution is skewed up to $a_{\mathrm{b}}\sim
650$ AU, comparable to the expected distribution for equal-mass ratio systems
(see center panel of Fig. 5 and Sect. 5.1).

Because a more massive companion increases the magnitude of the binary’s
torque, spin-orbit misalignments occur more often in circum-secondary than
circum-primary systems. Few circum-primary systems experience obliquity
excitations, and in general only within binary separations of $\lesssim 500$
AU. On the other hand, circum-secondary systems almost always have misaligned
stellar obliquities, even at larger binary separations.

##  6 The Mass Dependence of Alignment in Observed Systems

Figure 8: Stellar mass ratio vs binary separation (left panel) and mass ratio
vs sky-projected inclination (center and right panel) for the sample of
planet-hosting binaries selected and curated by Christian et al. (2024). We
delineate the sample into two subsets lying inside and outside a binary
separation of $800$ AU, as this is approximately where Christian et al. (2022)
find evidence for preferential alignment. Figure 9: Observed distributions of
sky-projected inclinations for circum-primary planets (left panel), near-equal
mass ratio systems (center panel), and circum-secondary planets (right panel).
The data are split into two semi-major axis bins inside or outside 800 AU. The
control sample is shown in dotted lines. The fact that the control sample
deviates from isotropy is indicative of systematic biases as discussed in
Christian et al. (2022). We also show the $p$-value resulting from a KS-test
for each population.

The damping via dissipative precession depends on stellar mass ratio, as
discussed in Sects. 5.1 and 5.2. It is thus of immediate interest to
investigate a potential correlation between stellar mass ratio and orbital
orientations in observed exoplanetary populations.

###  6.1 Sample selection and data curation

Our sample is drawn from the intersection of systems observed with one or more
transiting exoplanet candidates in the Gaia DR2-based TESS Input Catalog
(TIC); two resolved stars in Gaia DR3; and a sufficiently robust astrometric
solution/low probability of chance alignment as quantified by El-Badry et al.
(2021). After removal of false positives, the sample includes 398 candidate
exoplanet-hosting binaries. We also leverage the control sample of binary
systems without observed transiting exoplanets from Christian et al. (2024) in
our analysis. This allows us to filter for selection effects. Stellar masses
for both populations were obtained using a combination of empirical mass-
luminosity relations and stellar isochrone fitting. We refer to Christian et
al. (2024) and Christian et al. (2022) for further details on the data
selection and curation, as well as the control sample.

For all of the binaries in the sample, the separation is less than 21
arcseconds on the sky, which is the size of a single TESS pixel. It is
therefore difficult in general to determine which of the two stars the
transiting exoplanet is orbiting. While in some cases the TESS data can be
used to identify the true host star by fitting the telescope’s pixel response
function to the in/out-of-transit difference image, in many other cases
(especially for the closest binaries that are most likely to be aligned) we
cannot confidently identify which star is the host unless there is ground-
based follow-up data.

In order to proceed, we assume that all planets in our sample orbit the host
stars indicated in the TESS Objects of Interest catalog (Guerrero et al.,
2021). Effectively, this means that we assume planets all orbit the more
massive (brighter) star unless there is compelling evidence to the contrary
(either from TESS centroid analysis or ground-based follow-up observations).
This is generally a good assumption because of the biases of the transit
method, which is more sensitive to deeper transit signals that are more easily
observable from the brighter star of a blended signal (Ciardi et al., 2015).

Figure 8 shows stellar mass ratio vs binary separation for the sample (left
panel), as well as mass ratio vs measured inclination (center and right
panel). There are 15 systems with a mass ratio of $M_{*}/M_{\mathrm{b}}\leq
0.8$, only four of which have semimajor axes less than 800 AU. This is too few
for any meaningful statistical statements. On the other hand, there is a large
number of systems that we can broadly assign to the equal mass ratio and
circum-primary categories, defined here as $0.8<M_{*}/M_{\mathrm{b}}\leq 2$
and $M_{*}/M_{\mathrm{b}}>2$, respectively. The data thus present an ideal
proving ground to test the significance of the difference in the inclination
distribution between equal-mass ratio and circum-primary systems (see Fig. 5).

###  6.2 Orientation distribution as a function of stellar mass ratio

Figure 9 shows the $\cos(i_{\mathrm{b}})$ distribution of the sample,
separated by stellar mass ratio and binary separation. We overplot equivalent
distributions for the control sample, which includes no transiting planets,
yet likewise deviates from isotropy due to systematic biases (see Christian et
al., 2022, for a discussion). It is thus more appropriate to compare the
planet-hosting sample to the control sample rather than the line of isotropy
which is drawn in Fig. 9 for reference.

Equal mass ratio systems (center panel of Fig. 9) exhibit a trend consistent
with that published in Christian et al. (2022) and Christian et al. (2024):
exoplanet-hosting binaries with semi-major axes greater than $\sim$ 800 AU are
consistent with the (almost isotropic) control sample, whereas those with
semi-major axes less than $\sim$ 800 AU visibly display preferential
alignment. This is confirmed by a two-sample Kolmogorov-Smirnov (KS) test,
where the resulting $p$-value characterizes the probability that the two
samples originate from the same underlying distribution. The null hypothesis
typically is rejected if $p<5\%$. Indeed, the two sample KS-test yields
$p$-values of $0.01\%$ and $80.4\%$ for the equal mass ratio populations
inside and outside 800 AU respectively, when comparing to the associated
control sample population.

The high-mass ratio systems are distributed – at least by eye – more
isotropically. Indeed, a KS test with the corresponding control sample yields
$p$-values of $29.4\%$ and $58.4\%$ for the high mass ratio populations inside
and outside 800 AU, respectively. The KS test therefore cannot reject the null
hypothesis that exoplanet-hosting high-mass ratio binaries are oriented in a
similar fashion to the binaries in the control sample.

To conclude, while the observational data are relatively sparse, they are
broadly consistent with the theoretical predictions made in Sect. 5.2 and Fig.
5.

##  7 Additional alignment pathways

In our analysis of exoplanet-hosting binary systems, we observed
characteristics consistent with alignment through viscous dissipation.
However, our findings also indicate that viscous dissipation alone may not
fully account for the current orientation distribution of these systems (see
Sect. 5.1). While present-day data is ambiguous, if polar systems are indeed
rare, this phenomenon could not be adequately reproduced by dissipative
precession – at least as treated in the models presented here.

In this section, we therefore discuss two other possibilities – namely,
preferential primordial alignment and removal of highly inclined systems
through the Kozai-Lidov mechanism – alongside their prospects to, jointly with
dissipative precession, sculpt the orbital configurations of planet-hosting
binaries.

###  7.1 Primordial alignment & magnetic torques

Our models assumed an initial configuration in which the spin and disk angular
momentum vectors were aligned but randomly oriented relative to the binary
angular momentum vector, leading to a probability density function for the
invariant disk – binary inclination $\theta_{\mathrm{db}}$ given by
$\sin|\theta_{\mathrm{db}}|$. As discussed in Sect. 4.1, this choice aims to
reflect turbulent fragmentation, where the binary components are assumed to
form independently and with random orientations, as the mode of star formation
examined within this work.

Numerical investigations, however, suggest that, even for wide binary
formation via turbulent fragmentation, a dearth of polar systems is expected
(e.g., Bate et al., 2010; Bate, 2018). Physically, this is a consequence of
the conservation of the angular momentum of the contracting cloud core, which,
even in the presence of strong turbulence, will lead to a non-isotropic
orientation of the angular momentum vectors of individual stars.

In addition, accretion of ambient gas onto the binary components acts to align
the binary inclination (Bate, 2018). Torques induced by gas accretion can also
hinder excitation of spin-orbit misalignment (Lai, 2014), thus affecting
expected obliquity distributions beyond what is discussed in our paper. The
obliquity distribution can, furthermore, be impacted by magnetic torques. Lai
et al. (2011); Lai (2014) show that magnetic torques contribute to precession
rates and thus to the excitation of $\theta_{\mathrm{sd}}$. On the other hand,
magnetic torques also can act to (re-)align stellar spin with the disk,
specifically for low mass stars (Spalding & Batygin, 2015). As magnetic
effects generally do not affect the mutual disk-binary inclination – the focus
of this paper – we do not include them in this study.

###  7.2 Kozai-Lidov oscillations in polar configurations

Kozai-Lidov oscillations are suppressed as long as the disk self-gravity is
strong enough to force eccentric orbits in apsidal precession (see Sect 3.3).
On the other hand, once the protoplanetary disk evaporates, and provided the
binary inclination is sufficiently large, the Kozai-Lidov mechanism can freely
operate, in the process affecting the orbital configurations of the planets
within the system. Specifically, this can lead to a distribution where polar
configurations are removed (Naoz et al., 2012). Planet-planet interactions
become more frequent in systems undergoing Kozai-Lidov cycles (Malmberg et
al., 2007), which, for highly inclined systems, may lead to scattering and
thus the removal of such systems from the observed orientation distributions.

The Kozai-Lidov mechanism may, therefore, reconcile the inclination
distributions produced by this paper’s model (e.g., Fig. 5) with the observed
distributions in Fig. 9 or in Christian et al. (2022); Rice et al. (2024),
where polar systems are uncommon.

##  8 Conclusions

In this paper, we investigated the prospects for viscous dissipation –
occurring during precession of the protoplanetary disk driven by an inclined
binary companion – to sculpt resulting orbital configurations. Physically, the
disk’s mass loss is associated with a resonant crossing of precession
frequencies which excites the system’s spin-orbit angle (Sect. 2). At the same
time, the orbit-orbit angle is damped as the disk dissipates energy in an
effort to resist twisting by the torque from the binary companion (Sect. 3).
The sky-projected inclination and obliquity distributions predicted by this
model are broadly consistent with observed trends (Sects. 4 and 5).
Specifically, if dissipative precession can operate effectively, it may offer
a joint explanation for:

  * •

The preferential population-level alignment of binaries that host transiting
exoplanets perpendicular to the sky plane, and the decrease in alignment for
increasing binary separations (Dupuy et al., 2022; Christian et al., 2022).

  * •

The novel finding that this preferential alignment appears to be more
pronounced for binaries with close-to equal stellar mass ratios compared to
circum-primary planets within high mass ratio binaries, where it may not exist
at all (Sect. 6).

  * •

The stellar obliquity distribution, which consists of mostly aligned systems,
with a second population of significantly misaligned or isotropically
distributed systems (Albrecht et al., 2021; Dong & Foreman-Mackey, 2023;
Siegel et al., 2023; Rice et al., 2024).

An important caveat is that the damping rate due to viscous dissipation during
binary-driven disk precession weakens with increasing binary separation
$a_{\mathrm{b}}$ and decreasing disk extent $r_{\mathrm{out}}$, viscosity
$\alpha$ and binary mass ratio $M_{\mathrm{b}}/M_{*}$ via Eq. (26). Observed
protoplanetary disks show a large diversity in these properties (see Sect.
3.5), and consequently, dissipative precession is expected to only be
efficient for the sub-set of systems adequately described with parameters
around those chosen by our fiducial model ($r_{\mathrm{out}}\gtrsim 100$ AU,
$\alpha\gtrsim 0.05,M_{\mathrm{b}}/M_{*}\gtrsim 1$). Additional constraints on
disk properties, specifically their radial extent and their viscosity, will
help to pinpoint the role of alignment via dissipative precession.

Our model presupposes that star formation results in a spherically isotropic
distribution of initial inclinations. Such a distribution would leave behind a
population of polar systems seemingly absent in the observational data.
However, this discrepancy can be mitigated by viscous dissipation during
binary-driven recession juxtaposed with even small degrees of primordial
alignments, removal of polar systems via Kozai-Lidov oscillations, or both
(Sect. 7). While the presented model can reproduce the broad, qualitative
trends under favorable conditions, we suspect that a full sophisticated
population synthesis, which jointly models these effects and covers the full
parameter space of appropriate disk properties, is required to match the
observed, joint - obliquity - binary inclination - distribution in Rice et al.
(2024). As the observational sample of obliquity measurements in exoplanet-
hosting binaries grows, so will the appeal of such an undertaking.

Finally, our work demonstrates that the distributions of planet orbit
orientations and stellar obliquities in binaries may be dependent on stellar
mass ratio, thus highlighting that planet-hosting binaries are a unique window
into the planet formation process. In this specific context, additional
demographic constraints on planets in high mass-ratio binary systems, in
particular, difficult-to-identify circum-secondary planets, may be
exceptionally decisive in diagnosing the predominant mechanism that drives
exoplanets towards alignment with their stellar binary.

##  9 Acknowledgements

We thank the referee for insightful comments. K.G. thanks Jeremy Smallwood,
Mark Dodici, Yubo Su, Jiayin Dong, Callum Fairbairn and Tiger Lu for
insightful discussions. M.R. acknowledges support from Heising-Simons
Foundation Grants #2023-4478 and #2023-4655, NASA grant GR122985/AWD0011072,
and Oracle for Research grant No. CPQ-3033929. JJZ was supported by a 51
Pegasi b Heising-Simons Fellowship.

.

##  Appendix A Calculation of the relative longitude $\Delta\Omega$

Consider the set of vectors and angles depicted in Fig. 1. We denote the
projections of $\hat{\bm{s}}$ and $\hat{\bm{l}}_{\mathrm{b}}$ onto the $x-z$
plane as $\hat{\bm{s}}_{xz}$ and $\hat{\bm{l}}_{\mathrm{b},xz}$ respectively,
and the projections onto the $y$-axis as $\hat{\bm{s}}_{y}$ and
$\hat{\bm{l}}_{\mathrm{b},y}$ respectively.

$\hat{\bm{s}}$ and $\hat{\bm{l}}_{\mathrm{b}}$ create a triangle, where the
remaining side $\overline{SB}$ is opposite to the angle
$\theta_{\mathrm{sb}}$. Since $\hat{\bm{s}}$ and $\hat{\bm{l}}_{\mathrm{b}}$
are unit vectors, the law of cosines implies

$\displaystyle\overline{SB}^{2}=2(1-\cos\theta_{\mathrm{sb}}).$ (A1)

We continue by considering the plane $P$ parallel to the $y$-axis that
contains $\overline{SB}$. At $y=0$, the plane also contains the remaining side
$\overline{Y}$ of the triangle constructed by $\hat{\bm{s}}_{xz}$ and
$\hat{\bm{l}}_{\mathrm{b},xz}$. $\overline{Y}$ is opposite to the relative
longitude $\Delta\Omega$, so we apply the law of cosines again as

$\displaystyle\cos\Delta\Omega=\frac{|\hat{\bm{s}}_{xz}|^{2}+|\hat{\bm{l}}_{\mathrm{b},xz}|^{2}-\overline{Y}^{2}}{2|\hat{\bm{s}}_{xz}||\hat{\bm{l}}_{\mathrm{b},xz}|},$
(A2)

where $|\hat{\bm{s}}_{xz}|=\sin\theta_{\mathrm{sd}}$ and
$|\hat{\bm{l}}_{\mathrm{b},xz}|=\sin\theta_{\mathrm{db}}$ since $\hat{\bm{s}}$
and $\hat{\bm{l}}_{\mathrm{b}}$ are unit vectors, and
$\hat{\bm{l}}_{\mathrm{d}}$ is aligned with the $y$-axis. We consider a right-
angled triangle in $P$, that has hypotenuse $\overline{SB}$ and the legs
$\overline{Y}$ and $\overline{D}$, such that
$\overline{SB}^{2}=\overline{Y}^{2}+\overline{D}^{2}$. Here, $\overline{D}$ is
parallel to $\hat{\bm{l}}_{\mathrm{b}}$ and the $y$-axis, and can be defined
as $\overline{D}=|\hat{\bm{s}}_{y}|-|\hat{\bm{l}}_{\mathrm{b},y}|$. These
projections have magnitudes of $|\hat{\bm{s}}_{y}|=\cos\theta_{\mathrm{sd}}$
and $|\hat{\bm{l}}_{\mathrm{b},y}|=\cos\theta_{\mathrm{db}}$.

We can therefore write

$\displaystyle\overline{Y}^{2}=\overline{SB}^{2}-\overline{D}^{2}=\overline{SB}^{2}-(|\hat{\bm{s}}_{y}|-|\hat{\bm{l}}_{\mathrm{b},y}|)^{2}=2(1-\cos\theta_{\mathrm{sb}})-(\cos\theta_{\mathrm{sd}}-\cos\theta_{\mathrm{db}})^{2},$
(A3)

and thus

$\displaystyle\cos\Delta\Omega$
$\displaystyle=\frac{\sin^{2}\theta_{\mathrm{db}}+\sin^{2}\theta_{\mathrm{sd}}+(\cos\theta_{\mathrm{sd}}-\cos\theta_{\mathrm{db}})^{2}-2(1-\cos\theta_{\mathrm{sb}})}{2\sin\theta_{\mathrm{db}}\sin\theta_{\mathrm{sb}}},$
(A4)

which simplifies to Eq. (35).

## References

  * Albrecht et al. (2021) Albrecht, S. H., Marcussen, M. L., Winn, J. N., Dawson, R. I., & Knudstrup, E. 2021, The Astrophysical Journal Letters, 916, L1
  * Andersen et al. (2019) Andersen, B. C., Stephens, I. W., Dunham, M. M., et al. 2019, ApJ, 873, 54, doi: 10.3847/1538-4357/ab05c7
  * Anderson & Lai (2018) Anderson, K. R., & Lai, D. 2018, MNRAS, 480, 1402, doi: 10.1093/mnras/sty1937
  * Anderson et al. (2016) Anderson, K. R., Storch, N. I., & Lai, D. 2016, MNRAS, 456, 3671, doi: 10.1093/mnras/stv2906
  * Andrews & Williams (2005) Andrews, S. M., & Williams, J. P. 2005, ApJ, 631, 1134, doi: 10.1086/432712
  * Ansdell et al. (2018) Ansdell, M., Williams, J. P., Trapman, L., et al. 2018, ApJ, 859, 21, doi: 10.3847/1538-4357/aab890
  * Bate (2018) Bate, M. R. 2018, MNRAS, 475, 5618, doi: 10.1093/mnras/sty169
  * Bate et al. (2010) Bate, M. R., Lodato, G., & Pringle, J. E. 2010, MNRAS, 401, 1505, doi: 10.1111/j.1365-2966.2009.15773.x
  * Batygin (2012) Batygin, K. 2012, Nature, 491, 418, doi: 10.1038/nature11560
  * Batygin & Adams (2013) Batygin, K., & Adams, F. C. 2013, ApJ, 778, 169, doi: 10.1088/0004-637X/778/2/169
  * Batygin et al. (2011) Batygin, K., Morbidelli, A., & Tsiganis, K. 2011, A&A, 533, A7, doi: 10.1051/0004-6361/201117193
  * Birnstiel et al. (2012) Birnstiel, T., Klahr, H., & Ercolano, B. 2012, A&A, 539, A148, doi: 10.1051/0004-6361/201118136
  * Birnstiel et al. (2018) Birnstiel, T., Dullemond, C. P., Zhu, Z., et al. 2018, ApJ, 869, L45, doi: 10.3847/2041-8213/aaf743
  * Bonnell & Bate (1994) Bonnell, I. A., & Bate, M. R. 1994, MNRAS, 269, L45, doi: 10.1093/mnras/269.1.L45
  * Booth & Clarke (2019) Booth, R. A., & Clarke, C. J. 2019, MNRAS, 483, 3718, doi: 10.1093/mnras/sty3340
  * Boyden & Eisner (2020) Boyden, R. D., & Eisner, J. A. 2020, ApJ, 894, 74, doi: 10.3847/1538-4357/ab86b7
  * Chachan & Lee (2023) Chachan, Y., & Lee, E. J. 2023, arXiv e-prints, arXiv:2305.00803, doi: 10.48550/arXiv.2305.00803
  * Chiang & Goldreich (1997) Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368, doi: 10.1086/304869
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102, doi: 10.3847/0004-637X/823/2/102
  * Christian et al. (2022) Christian, S., Vanderburg, A., Becker, J., et al. 2022, AJ, 163, 207, doi: 10.3847/1538-3881/ac517f
  * Christian et al. (2024) —. 2024, arXiv e-prints, arXiv:2405.10379.  https://arxiv.org/abs/2405.10379
  * Ciardi et al. (2015) Ciardi, D. R., Beichman, C. A., Horch, E. P., & Howell, S. B. 2015, ApJ, 805, 16, doi: 10.1088/0004-637X/805/1/16
  * Dong & Foreman-Mackey (2023) Dong, J., & Foreman-Mackey, D. 2023, AJ, 166, 112, doi: 10.3847/1538-3881/ace105
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8, doi: 10.3847/0067-0049/222/1/8
  * Doğan et al. (2015) Doğan, S., Nixon, C., King, A., & Price, D. J. 2015, MNRAS, 449, 1251, doi: 10.1093/mnras/stv347
  * Doǧan et al. (2018) Doǧan, S., Nixon, C. J., King, A. R., & Pringle, J. E. 2018, MNRAS, 476, 1519, doi: 10.1093/mnras/sty155
  * Dullemond (2000) Dullemond, C. P. 2000, A&A, 361, L17, doi: 10.48550/arXiv.astro-ph/0007399
  * Dullemond et al. (2022) Dullemond, C. P., Kimmig, C. N., & Zanazzi, J. J. 2022, MNRAS, 511, 2925, doi: 10.1093/mnras/stab2791
  * Dupuy et al. (2022) Dupuy, T. J., Kraus, A. L., Kratter, K. M., et al. 2022, MNRAS, 512, 648, doi: 10.1093/mnras/stac306
  * Dupuy et al. (2022) Dupuy, T. J., Kraus, A. L., Kratter, K. M., et al. 2022, Monthly Notices of the Royal Astronomical Society, 512, 648
  * El-Badry et al. (2021) El-Badry, K., Rix, H.-W., & Heintz, T. M. 2021, MNRAS, 506, 2269, doi: 10.1093/mnras/stab323
  * Ercolano et al. (2009) Ercolano, B., Clarke, C. J., & Drake, J. J. 2009, ApJ, 699, 1639, doi: 10.1088/0004-637X/699/2/1639
  * Esposito et al. (2014) Esposito, M., Covino, E., Mancini, L., et al. 2014, A&A, 564, L13, doi: 10.1051/0004-6361/201423735
  * Fabrycky & Winn (2009) Fabrycky, D. C., & Winn, J. N. 2009, ApJ, 696, 1230, doi: 10.1088/0004-637X/696/2/1230
  * Flock et al. (2017) Flock, M., Nelson, R. P., Turner, N. J., et al. 2017, ApJ, 850, 131, doi: 10.3847/1538-4357/aa943f
  * Foucart & Lai (2014) Foucart, F., & Lai, D. 2014, MNRAS, 445, 1731, doi: 10.1093/mnras/stu1869
  * Fu et al. (2015) Fu, W., Lubow, S. H., & Martin, R. G. 2015, ApJ, 807, 75, doi: 10.1088/0004-637X/807/1/75
  * Gaia Collaboration et al. (2022) Gaia Collaboration, Vallenari, A., Brown, A. G. A., et al. 2022, arXiv e-prints, arXiv:2208.00211.  https://arxiv.org/abs/2208.00211
  * Gammie (2001) Gammie, C. F. 2001, ApJ, 553, 174, doi: 10.1086/320631
  * Gammie et al. (2000) Gammie, C. F., Goodman, J., & Ogilvie, G. I. 2000, MNRAS, 318, 1005, doi: 10.1046/j.1365-8711.2000.03669.x
  * Gerbig et al. (2019) Gerbig, K., Lenz, C. T., & Klahr, H. 2019, A&A, 629, A116, doi: 10.1051/0004-6361/201935278
  * Guerrero et al. (2021) Guerrero, N. M., Seager, S., Huang, C. X., et al. 2021, ApJS, 254, 39, doi: 10.3847/1538-4365/abefe1
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357, doi: 10.1038/s41586-020-2649-2
  * Hartmann et al. (1998) Hartmann, L., Calvet, N., Gullbring, E., & D’Alessio, P. 1998, ApJ, 495, 385, doi: 10.1086/305277
  * Hendler et al. (2020) Hendler, N., Pascucci, I., Pinilla, P., et al. 2020, ApJ, 895, 126, doi: 10.3847/1538-4357/ab70ba
  * 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
  * Hughes & Armitage (2012) Hughes, A. L. H., & Armitage, P. J. 2012, MNRAS, 423, 389, doi: 10.1111/j.1365-2966.2012.20892.x
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Komaki et al. (2021) Komaki, A., Nakatani, R., & Yoshida, N. 2021, ApJ, 910, 51, doi: 10.3847/1538-4357/abe2af
  * Kozai (1962) Kozai, Y. 1962, AJ, 67, 591, doi: 10.1086/108790
  * Kratter & Murray-Clay (2011) Kratter, K. M., & Murray-Clay, R. A. 2011, ApJ, 740, 1, doi: 10.1088/0004-637X/740/1/1
  * Krumholz et al. (2007) Krumholz, M. R., Klein, R. I., & McKee, C. F. 2007, ApJ, 656, 959, doi: 10.1086/510664
  * Kuruwita & Haugbølle (2023) Kuruwita, R. L., & Haugbølle, T. 2023, A&A, 674, A196, doi: 10.1051/0004-6361/202244882
  * 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
  * Larwood & Papaloizou (1997) Larwood, J. D., & Papaloizou, J. C. B. 1997, MNRAS, 285, 288, doi: 10.1093/mnras/285.2.288
  * Lee et al. (2017) Lee, J.-E., Lee, S., Dunham, M. M., et al. 2017, Nature Astronomy, 1, 0172, doi: 10.1038/s41550-017-0172
  * Lester et al. (2023) Lester, K. V., Howell, S. B., Matson, R. A., et al. 2023, AJ, 166, 166, doi: 10.3847/1538-3881/acf563
  * Li (2021) Li, G. 2021, ApJ, 915, L2, doi: 10.3847/2041-8213/ac0620
  * Lidov (1962) Lidov, M. L. 1962, Planet. Space Sci., 9, 719, doi: 10.1016/0032-0633(62)90129-0
  * Liu et al. (2022) Liu, Y., Linz, H., Fang, M., et al. 2022, A&A, 668, A175, doi: 10.1051/0004-6361/202244505
  * Lodato & Price (2010) Lodato, G., & Price, D. J. 2010, MNRAS, 405, 1212, doi: 10.1111/j.1365-2966.2010.16526.x
  * Lubow & Ogilvie (2000) Lubow, S. H., & Ogilvie, G. I. 2000, ApJ, 538, 326, doi: 10.1086/309101
  * Lubow & Ogilvie (2017) —. 2017, MNRAS, 469, 4292, doi: 10.1093/mnras/stx990
  * Malmberg et al. (2007) Malmberg, D., Davies, M. B., & Chambers, J. E. 2007, MNRAS, 377, L1, doi: 10.1111/j.1745-3933.2007.00291.x
  * Manara et al. (2023) Manara, C. F., Ansdell, M., Rosotti, G. P., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 539, doi: 10.48550/arXiv.2203.09930
  * Martin et al. (2014) Martin, R. G., Nixon, C., Lubow, S. H., et al. 2014, ApJ, 792, L33, doi: 10.1088/2041-8205/792/2/L33
  * Millholland & Batygin (2019) Millholland, S., & Batygin, K. 2019, ApJ, 876, 119, doi: 10.3847/1538-4357/ab19be
  * Miotello et al. (2023) Miotello, A., Kamp, I., Birnstiel, T., Cleeves, L. C., & Kataoka, A. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 501, doi: 10.48550/arXiv.2203.09818
  * Miranda & Lai (2015) Miranda, R., & Lai, D. 2015, MNRAS, 452, 2396, doi: 10.1093/mnras/stv1450
  * Najita & Bergin (2018) Najita, J. R., & Bergin, E. A. 2018, ApJ, 864, 168, doi: 10.3847/1538-4357/aad80c
  * Naoz (2016) Naoz, S. 2016, ARA&A, 54, 441, doi: 10.1146/annurev-astro-081915-023315
  * Naoz et al. (2012) Naoz, S., Farr, W. M., & Rasio, F. A. 2012, ApJ, 754, L36, doi: 10.1088/2041-8205/754/2/L36
  * Nixon & King (2012) Nixon, C. J., & King, A. R. 2012, MNRAS, 421, 1201, doi: 10.1111/j.1365-2966.2011.20377.x
  * Offner et al. (2023) Offner, S. S. R., Moe, M., Kratter, K. M., et al. 2023, in Astronomical Society of the Pacific Conference Series, Vol. 534, Protostars and Planets VII, ed. S. Inutsuka, Y. Aikawa, T. Muto, K. Tomida, & M. Tamura, 275, doi: 10.48550/arXiv.2203.10066
  * Ogilvie (1999) Ogilvie, G. I. 1999, MNRAS, 304, 557, doi: 10.1046/j.1365-8711.1999.02340.x
  * Ogilvie & Latter (2013) Ogilvie, G. I., & Latter, H. N. 2013, MNRAS, 433, 2420, doi: 10.1093/mnras/stt917
  * Owen et al. (2012) Owen, J. E., Clarke, C. J., & Ercolano, B. 2012, MNRAS, 422, 1880, doi: 10.1111/j.1365-2966.2011.20337.x
  * Owen et al. (2011) Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13, doi: 10.1111/j.1365-2966.2010.17818.x
  * Padoan & Nordlund (2002) Padoan, P., & Nordlund, Å. 2002, ApJ, 576, 870, doi: 10.1086/341790
  * Papaloizou & Lin (1995) Papaloizou, J. C. B., & Lin, D. N. C. 1995, ApJ, 438, 841, doi: 10.1086/175127
  * Papaloizou & Pringle (1983) Papaloizou, J. C. B., & Pringle, J. E. 1983, MNRAS, 202, 1181, doi: 10.1093/mnras/202.4.1181
  * Powell et al. (2019) Powell, D., Murray-Clay, R., Pérez, L. M., Schlichting, H. E., & Rosenthal, M. 2019, ApJ, 878, 116, doi: 10.3847/1538-4357/ab20ce
  * Rafikov (2017) Rafikov, R. R. 2017, ApJ, 837, 163, doi: 10.3847/1538-4357/aa6249
  * Rice et al. (2024) Rice, M., Gerbig, K., & Vanderburg, A. 2024, AJ, 167, 126, doi: 10.3847/1538-3881/ad1bed
  * Rice et al. (2023) Rice, M., Wang, S., Gerbig, K., et al. 2023, AJ, 165, 65, doi: 10.3847/1538-3881/aca88e
  * Rosotti (2023) Rosotti, G. P. 2023, New A Rev., 96, 101674, doi: 10.1016/j.newar.2023.101674
  * 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 et al. (2013) Sanchis-Ojeda, R., Winn, J. N., Marcy, G. W., et al. 2013, ApJ, 775, 54, doi: 10.1088/0004-637X/775/1/54
  * Shakura & Sunyaev (1973) Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337
  * Siegel et al. (2023) Siegel, J. C., Winn, J. N., & Albrecht, S. H. 2023, ApJ, 950, L2, doi: 10.3847/2041-8213/acd62f
  * Smallwood et al. (2023) Smallwood, J. L., Martin, R. G., & Lubow, S. H. 2023, MNRAS, 520, 2952, doi: 10.1093/mnras/stad338
  * Spalding & Batygin (2015) Spalding, C., & Batygin, K. 2015, ApJ, 811, 82
  * Spalding et al. (2014) Spalding, C., Batygin, K., & Adams, F. C. 2014, ApJ, 797, L29, doi: 10.1088/2041-8205/797/2/L29
  * Tobin et al. (2018) Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2018, ApJ, 867, 43, doi: 10.3847/1538-4357/aae1f7
  * Tokovinin (2017) Tokovinin, A. 2017, MNRAS, 468, 3461, doi: 10.1093/mnras/stx707
  * van der Velden (2020) van der Velden, E. 2020, The Journal of Open Source Software, 5, 2004, doi: 10.21105/joss.02004
  * Vick et al. (2023) Vick, M., Su, Y., & Lai, D. 2023, ApJ, 943, L13, doi: 10.3847/2041-8213/acaea6
  * Villenave et al. (2020) Villenave, M., Ménard, F., Dent, W. R. F., et al. 2020, A&A, 642, A164, doi: 10.1051/0004-6361/202038087
  * von Zeipel (1910) von Zeipel, H. 1910, Astronomische Nachrichten, 183, 345, doi: 10.1002/asna.19091832202
  * Williams & Cieza (2011) Williams, J. P., & Cieza, L. A. 2011, ARA&A, 49, 67, doi: 10.1146/annurev-astro-081710-102548
  * Winn et al. (2010) Winn, J. N., Johnson, J. A., Howard, A. W., et al. 2010, ApJ, 723, L223, doi: 10.1088/2041-8205/723/2/L223
  * Wright (1987) Wright, E. L. 1987, ApJ, 320, 818, doi: 10.1086/165597
  * Wu & Lithwick (2021) Wu, Y., & Lithwick, Y. 2021, ApJ, 923, 123, doi: 10.3847/1538-4357/ac2b9c
  * Wu & Murray (2003) Wu, Y., & Murray, N. 2003, ApJ, 589, 605, doi: 10.1086/374598
  * Xiang-Gruess & Papaloizou (2014) Xiang-Gruess, M., & Papaloizou, J. C. B. 2014, MNRAS, 440, 1179, doi: 10.1093/mnras/stu308
  * Zanazzi & Lai (2017) Zanazzi, J. J., & Lai, D. 2017, MNRAS, 467, 1957, doi: 10.1093/mnras/stx208
  * Zanazzi & Lai (2018a) —. 2018a, MNRAS, 477, 5207, doi: 10.1093/mnras/sty951
  * Zanazzi & Lai (2018b) —. 2018b, MNRAS, 478, 835, doi: 10.1093/mnras/sty1075
  * Zhang et al. (2023) Zhang, J., Weiss, L. M., Huber, D., et al. 2023, arXiv e-prints, arXiv:2310.03299, doi: 10.48550/arXiv.2310.03299
  * Zhang et al. (2018) Zhang, Y., Li, Q., Xie, J.-W., et al. 2018, ApJ, 861, 116, doi: 10.3847/1538-4357/aac6c3