# Accelerating the Hypergeometric Function with the Beta Integral to Derive
New Infinite Series for $\pi$ and Values of the Gamma Function

Cetin Hakimoglu-Brown
<EMAIL_ADDRESS>
###### Abstract

The beta integral is applied to accelerate the hypergeometric function
$2F1\left\\{1,B;C;w\right\\}$ to derive new infinite series for constants such
as $\pi$ and values of the gamma function. A compendium of new infinite series
is given. Ramanujan-like formulas for pi are also derived based on elementary
inverse trigonometric functions, including a formula with rational values that
adds 2.5 digits per terms, which makes the series much more compact than
similar formulas in the existing literature.

Keywords: infinite series, hypergeometric function, $\pi$ formula, gamma
function, universal constants, catalan’s constant, convergence improvement,
binomial sums, beta function, pochhammer

##  1 Introduction

At around the ’90s and early 2000s, David & Peter Borwein et al. [2] Travis
Sherman [3], K. A. Driver et al. [13], Batir [14] Christian Krattenthaler et
al. [4], and Boris Gour’evitch & Jes’us Guillera Goyanes [1], Bellard [6]
derived various binomial-like infinite series for constants such as $\pi$, by
using the beta integral as applied to inverse trigonometric integrals. This
paper continues on such earlier work, and extends it to the Euler integral of
the hypergeometric function, in order to make the the algorithm more
comprehensive and generalized.

Examples of formulas that can be derived though the algorithm include the so-
called ’324 formula’ for the exceptionally fast rate of convergence:

$\pi=\frac{\sqrt{3}}{60}\sum_{n=0}^{\infty}\frac{(2n)!(130n+109)}{\left(\frac{7}{6}\right)_{n}\left(\frac{11}{6}\right)_{n}(-1296)^{n}}$
(1.1)

And

$\frac{7\cdot 3^{1/2}(\Gamma(1/3))^{3}}{\pi\cdot
2^{1/3}}=\sum_{n=0}^{\infty}\frac{\left(\frac{2}{3}\right)_{n}\left(\frac{1}{2}\right)_{n}(102n+59)}{\left(\frac{13}{12}\right)_{n}\left(\frac{19}{12}\right)_{n}(-288)^{n}}$
(1.2)

These formulas were derived with the following integrals which we shall prove
in this paper:

$\frac{z\cdot\Gamma(a+b+2)}{\Gamma(a+1)\Gamma(b+1)}\int_{0}^{1}\frac{x^{a}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx=\sum_{n=0}^{\infty}\frac{(a+1)_{kn}(b+1)_{sn}}{(a+b+2)_{(k+s)n}(z)^{n}}$
(1.3)

We’ll also prove the more comprehensive integral-polynomial relation. This is
the so called ’machinery’ behind these infinite series:

$\displaystyle\frac{\Gamma(a+1)\Gamma(b+1)}{Z\cdot\Gamma(a+b+2)}\sum_{n=0}^{\infty}\left(\frac{(a+1)_{kn}(b+1)_{sn}}{(a+b+2)_{(k+s)n}(z)^{n}}(w)\right)=$
(1.4)
$\displaystyle\int_{0}^{1}\frac{Q(x)x^{a}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx=\int_{0}^{1}\frac{x^{a}(1-x)^{b}}{P(x)}dx=c$
$\displaystyle
w=\left\\{\begin{array}[]{c}a_{p}\prod_{g=1}^{p}\left(\frac{a+g+kn}{a+b+g+1+(k+s)n}\right)+a_{p-1}\prod_{g=1}^{p-1}\left(\frac{a+g+kn}{a+b+g+1+(k+s)n}\right)+\\\
\ldots+a_{1}\left(\frac{a+1+kn}{a+b+2+(k+s)n}\right)+a_{0}\end{array}\right\\}$

This paper was originally composed in 2009 and published in its current form
to SSRN in 2021 [8]. The final section includes new results. Later similar
results involving binomials since the publication of this paper include John
M. Campbell et al., 2019 [9], Wenchang Chu, 2011 [10], Kunle Adegok et al.,
2023 [11]. The ’324 formula’ and the exact implementation of the algorithm
used here was hitherto unknown, and independently re-derived by Zhi-Wei Sun,
2023 [12].

##  2 Derivation of Algorithm

We will begin by deriving formula (1.3) and apply it to derive (1.4).

Consider the beta function:

$\int_{0}^{1}x^{kn+a}(1-x)^{sn+b}dx=\frac{\Gamma(kn+a+1)\Gamma(sn+b+1)}{\Gamma((k+s)n+a+b+2)}$

Then divide both sides by $z^{n}$ and take the summation

$\sum_{n=0}^{\infty}\left(\int_{0}^{1}\frac{x^{kn+a}(1-x)^{sn+b}}{z^{n}}dx\right)=\sum_{n=0}^{\infty}\frac{\Gamma(kn+a+1)\Gamma(sn+b+1)}{\Gamma(a+b+2+n(s+k))\cdot
z^{n}}$

This is equal to:

$\sum_{n=0}^{\infty}\left(\int_{0}^{1}x^{a}(1-x)^{b}\left(\frac{x^{k}(1-x)^{s}}{z}\right)^{n}dx\right)=\sum_{n=0}^{\infty}\frac{\Gamma(kn+a+1)\Gamma(sn+b+1)}{\Gamma(a+b+2+n(s+k))\cdot
z^{n}}$ (2.0)

Summing an infinite geometric series:

$\sum_{n=0}^{\infty}\left(\int_{0}^{1}x^{a}(1-x)^{b}\left(\frac{x^{k}(1-x)^{s}}{z}\right)^{n}dx\right)=z\int_{0}^{1}\frac{x^{a}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx$
(2.1)

Using the identity: $(\Gamma(w))(w)_{n}=\Gamma(w+n)$ we derive:

$\frac{\Gamma(kn+a+1)\Gamma(sn+b+1)}{\Gamma((k+s)n+a+b+2)}=\frac{(a+1)_{kn}(b+1)_{sn}\Gamma(a+1)\Gamma(b+1)}{(a+b+2)_{(k+s)n}\Gamma(a+b+2)}$
(2.2)

Plugging equations (2.1) and (2.2) into (2.0) and re-arranging terms we get
equation (1.3)

To make the integral useful for computing constants begin with a so called
’seed integral’. Such an integral must be in the form (for linear $P(x)$ is
the Euler integral for the hypergeoemtric function):

$\int_{0}^{1}\frac{x^{a}(1-x)^{b}}{P(x)}dx=c\text{ where }P(x)\text{ is a
polynomial and }c\text{ is a constant }$

Then choose integer values of $k$ and $s$ such that the polynomial expansion
of $z-x^{k}(1-x)^{s}$ is divisible by $P(x)$. The value of $z$ is determined
by computing the quotient polynomial $Q(x)$ via polynomial long division such
that the following identity is obtained:

$\int_{0}^{1}\frac{x^{a}(1-x)^{b}}{P(x)}dx=\int_{0}^{1}\frac{Q(x)x^{a}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx$

$Q(x)$ can be generalized to: $Q(x)=a_{p}x^{p}+a_{p-1}x^{p-1}\ldots
a_{1}x+a_{0}$

Multiplying $\int_{0}^{1}\frac{x^{a}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx$ by $Q(x)$
and splitting the $a_{p}$ coefficients of $Q(x)$ yields:

$\left\\{\begin{array}[]{l}a_{p}\int_{0}^{1}\frac{x^{a+p}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx+a_{p-1}\int_{0}^{1}\frac{x^{a+p-1}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx\\\
\ldots
a_{1}\int_{0}^{1}\frac{x^{a+1}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx+a_{0}\int_{0}^{1}\frac{x^{a}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx\end{array}\right\\}=\int_{0}^{1}\frac{x^{a}(1-x)^{b}}{P(x)}dx$

Each of these $a_{p}\ldots a_{0}$ split terms can be plugged into equation
(1.3) :

$\left\\{\begin{array}[]{l}a_{p}\int_{0}^{1}\frac{x^{a+p}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx=\frac{a_{p}\cdot\Gamma(a+p+1)\Gamma(b+1)}{z\cdot\Gamma(a+b+p+2)}\sum_{n=0}^{\infty}\frac{(a+p+1)_{kn}(b+1)_{sn}}{(a+b+p+2)_{(k+s)n}(z)^{n}}\\\
a_{p-1}\int_{0}^{1}\frac{x^{a+p-1}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx=\frac{a_{p-1}\cdot\Gamma(a+p)\Gamma(b+1)}{z\cdot\Gamma(a+b+p+1)}\sum_{n=0}^{\infty}\frac{(a+p)_{kn}(b+1)_{sn}}{(a+b+p+1)_{(k+s)n}(z)^{n}}\\\
\text{ - }\\\
a_{1}\int_{0}^{1}\frac{x^{a+1}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx=\frac{a_{1}\cdot\Gamma(a+2)\Gamma(b+1)}{z\cdot\Gamma(a+b+3)}\sum_{n=0}^{\infty}\frac{(a+2)_{kn}(b+1)_{sn}}{(a+b+3)_{(k+s)n}(z)^{n}}\\\
a_{0}\int_{0}^{1}\frac{x^{a}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx=\frac{a_{0}\cdot\Gamma(a+1)\Gamma(b+1)}{z\cdot\Gamma(a+b+2)}\sum_{n=0}^{\infty}\frac{(a+1)_{kn}(b+1)_{sn}}{(a+b+2)_{(k+s)n}(z)^{n}}\end{array}\right\\}$
(2.3)

Also note the following pochhammer and gamma function identities:

$\displaystyle(a+b+p+2)_{(k+s)n}=\left(\prod_{g=1}^{p}\left(\frac{a+b+g+1+(k+s)n}{a+b+g+1}\right)\right)(a+b+2)_{(k+s)n}$
(2.4)
$\displaystyle(a+p+1)_{kn}=\left(\prod_{g=1}^{p}\left(\frac{a+g+kn}{a+g}\right)\right)(a+1)_{(kn)}$
$\displaystyle\Gamma(a+b+p+2)=\left(\prod_{g=1}^{p}(a+b+g+1)\right)\Gamma(a+b+2)$
$\displaystyle\Gamma(a+p+1)=\left(\prod_{g=1}^{p}(a+g)\right)\Gamma(a+1)$

Putting it all together using (2.3) and (2.4) gives:

$\displaystyle\frac{a_{p}\cdot\Gamma(a+p+1)\Gamma(b+1)}{z\cdot\Gamma(a+b+p+2)}\sum_{n=0}^{\infty}\frac{(a+p+1)_{kn}(b+1)_{sn}}{(a+b+p+2)_{(k+s)n}(z)^{n}}=$
$\displaystyle\frac{a_{p}\left(\prod_{g=1}^{p}(a+g)\right)\Gamma(a+1)\Gamma(b+1)}{z\cdot\left(\prod_{g=1}^{p}(a+b+g+1)\right)\Gamma(a+b+2)}\sum_{n=0}^{\infty}\left(\frac{\left(\prod_{g=1}^{p}\left(\frac{a+g+kn}{a+g}\right)\right)(a+1)_{kn}(b+1)_{sn}}{\left.\left(\prod_{g=1}^{p}\left(\frac{a+b+g+1+(k+s)n}{a+b+g+1}\right)\right)(a+b+2)_{(k+s)n}(z)^{n}\right)}\right)$

Also note how the $\prod_{g=1}^{p}(a+b+g+1)$ and $\prod_{g=1}^{p}(a+g)$ are
eliminated such that we simplify and obtain:

$\displaystyle\frac{a_{p}\Gamma(a+1)\Gamma(b+1)}{z\cdot\Gamma(a+b+2)}\sum_{n=0}^{\infty}\left(\frac{\left(\prod_{g=1}^{p}(a+g+kn)\right)(a+1)_{kn}(b+1)_{sn}}{\left(\prod_{g=1}^{p}(a+b+g+1+(k+s)n)\right)(a+b+2)_{(k+s)n}(z)^{n}}\right)=$
(2.5) $\displaystyle
a_{p}\int_{0}^{1}\frac{x^{a+p}(1-x)^{b}}{z-x^{k}(1-x)^{s}}dx$

Performing this procedure for $a_{p}\ldots a_{0}$, summing the $\mathrm{p}+1$
terms, and factoring out
$\frac{(a+1)_{kn}(b+1)_{sn}}{(a+b+2)_{(k+s)n}(z)^{n}}$ gives equation (1.4),
which completes the proof.

One final pochhammer identity that will be encountered frequently in this
paper is:

$(a)_{nk}=\left(\prod_{y=0}^{k-1}\left(\frac{a+y}{k}\right)_{n}\right)k^{kn}$

Another formula that will be used extensively later is the Euler integral:

$2F1\left\\{x_{1},x_{2};y_{1};z\right\\}=\frac{\Gamma\left(y_{1}\right)}{\Gamma\left(x_{2}\right)\Gamma\left(y_{1}-x_{2}\right)}\int_{0}^{1}\frac{x^{x_{2}-1}(1-x)^{y_{1}-x_{2}-1}}{(1-zx)^{x_{1}}}dx$
(2.6)

###  2.1 Deriving Formulas for Pi

We shall use (1.4) to prove (1.1). Consider a hypergeometric function for
arcsine converted into integral form via (2.6):

$2F1\left\\{1,\frac{1}{2};\frac{3}{2};\frac{w^{2}}{w^{2}-1}\right\\}=\left(\frac{\sin^{-1}(w)}{w}\right)\sqrt{\left(1-w^{2}\right)}=\frac{1}{2}\int_{0}^{1}\frac{x^{-1/2}}{1-\left(\frac{x\cdot
w^{2}}{w^{2}-1}\right)}dx$ (2.7)

Then let $\mathrm{w}=1/2$ so that
$\frac{\pi\sqrt{3}}{3}=\int_{0}^{1}\frac{x^{-1/2}}{1+\frac{x}{3}}dx$

Using (1.4) let $\mathrm{k}=1$ and $\mathrm{s}=2$. Then let
$\mathrm{a}=-1/2\mathrm{~{}b}=0$ and $P(x)=1+\frac{x}{3}$ and find
$\mathrm{Q}(\mathrm{x})$ such that:

$\int_{0}^{1}\frac{x^{-1/2}Q(x)}{z-x(1-x)^{2}}dx=\int_{0}^{1}\frac{x^{-1/2}}{1+\frac{x}{3}}dx$

Expand $z-x(1-x)^{2}$ and perform polynomial division
$\frac{x}{3}+1\overline{\smash{)}-x^{3}+2x^{2}-x+z}$ A value of
$\mathrm{z}=-48$ and $Q(x)=-3x^{2}+15x-48$ is obtained. Thus, from which (1.1)
follows, completing the proof:

$\int_{0}^{1}\frac{x^{-1/2}\left(x^{2}-5x+16\right)}{-48-x(1-x)^{2}}dx=\frac{-1}{3}\int_{0}^{1}\frac{x^{-1/2}}{1+\frac{x}{3}}dx=\frac{-\pi\sqrt{3}}{9}$
(2.8)

(2.7) can be generalized, in which (1.1) immediately follows for $w=1/2$ :

$\displaystyle\sum_{n=0}^{\infty}\left(\frac{4w^{6}}{27(w^{2}-1)}\right)^{n}\frac{(1)_{n}(1/2)_{n}\left(n(4w^{4}+6w^{2}-18)+2w^{4}+5w^{2}-15\right)}{(11/6)_{n}(7/6)_{n}}=$
(2.9)
$\displaystyle\left(-\frac{15\sin^{-1}(w)}{w}\right)\sqrt{\left(1-w^{2}\right)}$

This follows from:

$\frac{w^{2}}{1-w^{2}x}=\frac{x^{2}+x(1/w^{2}-1)+(1-w^{2})/w^{4}}{(w^{2}-1)/w^{6}-x^{2}(1-x)}$
(2.10)

Equation (1.1) converges at a rate of $\log(324)$ or 2.5 digits per term using
an arcsine formula for $1/2$, which to the best of my knowledge is faster than
any previously known formula of this type bi-3F2 series. The result is a
rather compact formula that converges much faster than similar formula in the
literature. Part of the reason it converges so fast is because two of the
pochammer symbols cancel out for (2.5) and fractions cancel out for (2.4)
owing to the factorization of 6 into two unique primes (e.g. $6n+3=2n+1$), so
instead of it being a 4F3 series, it’s a 3F2 series yet retains the
accelerated convergence rate.

Like similar binomial $\pi$ formulas, aesthetically it bears a resemblance to
the Ramanujan $\pi$ formulas, but is derived through elementary means instead
of elliptic integrals, and not the reciprocal of $\pi$. Like the Ramanujan
formulas, (2.9) can be split into the sum of two 3F2 hypergeometric series.

I conjecture 1.1 is the fastest $\pi$ formula that can be expressed in the
form $c\pi=\sum_{n=0}^{\infty}Q(n)(a+bn)z^{n}$ in which $c$ is algebraic and
$a,b,z$ are rational, and Q(n) is some array of pochammer symbols for a 3F2
hypergeoemtric series.

A faster formula can be obtained that adds five digits per term. Let
$s=4,k=2,a=-1/2$, and $\mathrm{b}=0.P(x)=\left(1+\frac{x}{3}\right)$ Then
expand $z-x^{2}(1-x)^{4}$ and performing division for $\mathrm{Q}(\mathrm{x})$
the following identity is obtained:

$\int_{0}^{1}\frac{x^{-1/2}\left(-x^{5}+7x^{4}-27x^{3}+85x^{2}-256x+768\right)}{2304-x^{2}(1-x)^{4}}dx=\frac{\pi\sqrt{3}}{9}$
(2.11)

After some labor, we have:

$\pi=\frac{\sqrt{3}}{6^{5}}\sum_{n=0}^{\infty}\left(\frac{[(4n)!]^{2}(6n)!}{(2n)!(12n)!(9)^{n+1}}\right)\left(\frac{127169}{12n+1}-\frac{1070}{12n+5}-\frac{131}{12n+7}+\frac{2}{12n+11}\right)$
(2.12)

To derive another pi formula let
$\mathrm{s}=2,\mathrm{k}=2,\mathrm{a}=-1/2,\mathrm{~{}b}=0$ and
$P(x)=1+\frac{x}{3}$ to obtain:

$\pi
2^{10}\sqrt{3}=\sum_{k=0}^{\infty}\frac{1}{9^{k}\left(\begin{array}[]{c}8k\\\
4k\end{array}\right)}\left(\frac{5717}{8k+1}-\frac{413}{8k+3}-\frac{45}{8k+5}+\frac{5}{8k+7}\right)$
(2.13)

This is similar to [4]

##  3 Multivariable Values of P(x) and Q(x) for Infinite Series

Binomial identities can be found using equations (1.3) and (1.4), but
extending it to multiple variables. In this section $\mathrm{Q}(\mathrm{x})$
becomes $\mathrm{Q}(\mathrm{x},\mathrm{w}),\mathrm{P}(\mathrm{x})$ becomes
$\mathrm{P}(\mathrm{x},\mathrm{w})$, and $\mathrm{z}$ becomes a function of
$\mathrm{w}$.

To begin, using equation (1.3) let
$\mathrm{a}=0,\mathrm{~{}b}=0,\mathrm{k}=1,\mathrm{~{}s}=1,\mathrm{z}=\mathrm{w}$
to obtain the identity:

$\int_{0}^{1}\frac{1}{w-x(1-x)}dx=\sum_{n=0}^{\infty}\frac{1}{w^{n+1}\left(\begin{array}[]{l}2n\\\
n\end{array}\right)(2n+1)}$ (3.1)

By evaluating the left hand integral we see that (3.1) is equal to [3]:

$\frac{4\tan^{-1}(1/\sqrt{4w-1})}{\sqrt{4w-1}}$ (3.2)

Do derive a faster series for (3.1) let
$\mathrm{k}=3,\mathrm{~{}s}=3,\mathrm{a}=0,\mathrm{~{}b}=0,\mathrm{z}=\mathrm{z}$,
and expand $\mathrm{z}-x^{3}(1-x)^{3}$

And let $\mathrm{P}(\mathrm{x},\mathrm{w})=x^{2}-x+w$. Performing polynomial
division, we need to find $\mathrm{Q}(\mathrm{x},\mathrm{w})$ and a new value
of $\mathrm{z}$ in terms of
$\mathrm{w}\cdot(x^{2}-x+w)\overline{\smash{)}(x^{6}-3x^{5}+3x^{4}-x^{3}+z)}$

We find: $Q(x,w)=\left(x^{4}-2x^{3}+(1-w)x^{2}+xw+w^{2}\right),z=w^{3}$
Therefore:

$\int_{0}^{1}\frac{\left(x^{4}-2x^{3}+(1-w)x^{2}+xw+w^{2}\right)}{w^{3}-x^{3}(1-x)^{3}}dx=\int_{0}^{1}\frac{1}{w-x(1-x)}dx$
(3.3)

With regard to equation (3.1) we can let $w=w^{-2}$,divide by $\mathrm{w}$,
and take the derivative in terms of w such that:

$\sum_{n=0}^{\infty}\frac{w^{2n+1}}{\left(\begin{array}[]{l}2n\\\
n\end{array}\right)(2n+1)}\frac{d}{dw}=\sum_{n=0}^{\infty}\frac{w^{2n}}{\left(\begin{array}[]{l}2n\\\
n\end{array}\right)}=\frac{w^{2}}{4-w^{2}}+\frac{4w\tan^{-1}\left(w/\sqrt{4-w^{2}}\right)}{\left(4-w^{2}\right)^{3/2}}$
(3.4)

An arbitrary number of integrals and derivatives of (3.4) can be taken in
terms of $\mathrm{w}$ to obtain various central binomial identities.

Applying formula (1.4) on (3.3), and after some labor obtain:

(The variable $r$ denotes the powers of $n$ that results from taking multiple
derivatives)

$\displaystyle 4\sum_{n=1}^{\infty}\frac{w^{n}}{\left(\begin{array}[]{l}2n\\\
n\end{array}\right)n^{r}}=$ (3.5) $\displaystyle 2w^{2}+\frac{2w^{2}}{3\cdot
2^{r}}+\sum_{n=1}^{\infty}\frac{w^{3n}}{\left(\begin{array}[]{l}6n\\\
3n\end{array}\right)}\left(\frac{w^{2}\left(9n^{2}+9n+2\right)}{(3n+2)^{r}(6n+3)(6n+1)}+\frac{2w(3n+1)}{(6n+1)(3n+1)^{r}}+\frac{4}{(3n)^{r}}\right)$

Let $\mathrm{r}=0$ and $w=w^{2}$ in (3.5). Then divide (3.5) by $\mathrm{w}$
and take multiple derivatives such that the denominators vanish, ensuring you
divide (3.5) by $\mathrm{w}$ each time before the derivative is taken. Then
finally let $\mathrm{w}=1$ to obtain the identity:

$3\sum_{n=1}^{\infty}\frac{63n^{2}-27n+4}{\left(\begin{array}[]{l}6n\\\
3n\end{array}\right)}=16\sum_{n=1}^{\infty}\frac{n^{2}}{\left(\begin{array}[]{l}2n\\\
n\end{array}\right)}-32\sum_{n=1}^{\infty}\frac{n}{\left(\begin{array}[]{l}2n\\\
n\end{array}\right)}+12\sum_{n=1}^{\infty}\frac{1}{\left(\begin{array}[]{l}2n\\\
n\end{array}\right)}$ (3.6)

Using (3.4) and taking multiple derivatives we arrive at:

$3\sum_{n=1}^{\infty}\frac{63n^{2}-27n+4}{\left(\begin{array}[]{l}6n\\\
3n\end{array}\right)}=\frac{40\pi\sqrt{3}}{81}+4$ (3.7)

A faster converting series can be obtained by letting
$\mathrm{k}=5,\mathrm{~{}s}=5$ such that

$\displaystyle\int_{0}^{1}\frac{Q(x,w)}{w^{5}-x^{5}(1-x)^{5}}dx=\int_{0}^{1}\frac{1}{P(x,w)}dx$
(3.8) $\displaystyle P(x,w)=w-x(1-x)$ $\displaystyle
Q(x,w)=\left(\begin{array}[]{l}x^{8}-4x^{7}+(6-w)x^{6}+(3w-4)x^{5}+\left(w^{2}-3w+1\right)x^{4}\\\
+\left(w-2w^{2}\right)x^{3}+\left(w^{2}-w^{3}\right)x^{2}+w^{3}x+w^{4}\end{array}\right)$

Using formula (1.4) and after some labor the identity is obtained:

$\displaystyle\sum_{n=0}^{\infty}\frac{213125n^{4}-278000n^{3}+139975n^{2}-26800n+1596}{\left(\begin{array}[]{l}10n\\\
5n\end{array}\right)}=$ (3.9) $\displaystyle
16\sum_{n=0}^{\infty}\frac{16n^{4}-128n^{3}+344n^{2}-352n+105}{\left(\begin{array}[]{l}2n\\\
n\end{array}\right)}=\frac{1120\pi\sqrt{3}}{81}+1728$

##  4 Hypergeometric Function Transformations for the Gamma Function

Infinite series for the gamma function such as (1.2) can be derived using
identities (1.3) and (1.4). In this section we’ll also prove a quadratic
hypergeometric transformation that is a byproduct of this method, as well as
other identities.

Begin with Kummer’s formula [5] and express it in integral form via (2.6):

$2F1\\{1,h;2-h;-1\\}=\frac{\Gamma(2-h)\Gamma(3/2)}{\Gamma(3/2-h)\Gamma(2)}=\frac{\Gamma(2-h)}{\Gamma(2-2h)\Gamma(h)}\int_{0}^{1}\frac{x^{h-1}(1-x)^{1-2h}}{1+x}dx$
(4.1)

Since $-2-x(1-x)=(x-2)(x+1)$ we can let
$\mathrm{k}=1,\mathrm{~{}s}=1,\mathrm{z}=-2$ and $Q(x)=x-2$ And after
isolating the integral in (4.1) we have the identity:

$\int_{0}^{1}\frac{x^{h-1}(1-x)^{1-2h}(x-2)}{-2-x(1-x)}dx=\int_{0}^{1}\frac{x^{h-1}(1-x)^{1-2h}}{x+1}dx=\frac{\Gamma(2-2h)\Gamma(h)\sqrt{\pi}}{2\Gamma(3/2-h)}$
(4.2)

The leftmost integral in (4.3) can be evaluated using (1.4). setting
$a_{0}=-2a_{1}=1$ $\mathrm{a}=\mathrm{h}-1\mathrm{~{}b}=1-2\mathrm{~{}h}$ :

$\frac{\Gamma(h)\Gamma(2-2h)}{-2\Gamma(2-h)}\sum_{n=0}^{\infty}\frac{(h)_{n}(2-2h)_{n}}{(2-h)_{2n}(-2)^{n}}\left(\frac{h+n}{2-h+2n}-2\right)=\frac{\sqrt{\pi}\Gamma(h)\Gamma(2-2h)}{2\Gamma(3/2-h)}$

Using (2.5), we obtain the series computing values of the gamma function:

$\frac{h(h+1)2^{2h-1}[\Gamma(h)]^{2}}{\Gamma(2h)}=\sum_{n=0}^{\infty}\frac{(1-h)_{n}(2h)_{n}(3n+3h+1)}{(-8)^{n}(h/2+1)_{n}(h/2+3/2)_{n}}$
(4.3)

To derive (1.2) let $\mathrm{h}=1/3$ in equation (4.3) and simplify:

$\frac{3^{1/2}2^{2/3}[\Gamma(1/3)]^{3}}{9\pi}=\sum_{n=0}^{\infty}\frac{(2/3)_{n}(2/3)_{n}(3n+2)}{(-8)^{n}(7/6)_{n}(5/3)_{n}}$

But $(2/3)_{n}(3n+2)=2(5/3)_{n}$ so we can write:

$\frac{3^{1/2}2^{2/3}[\Gamma(1/3)]^{3}}{18\pi}=\sum_{n=0}^{\infty}\frac{(2/3)_{n}}{(-8)^{n}(7/6)_{n}}=2F1\left\\{1,\frac{2}{3};\frac{7}{6};\frac{-1}{8}\right\\}$
(4.4)

(4.4)Is a quadratically transformed hypergeometric function evaluated at a
specific point. This to the best of my knowledge is a new derivation of this
$\Gamma(1/3)$ series, which does not involve either manipulation of the
hypergeometric function, transformation of the Euler integral, or evaluating
the hypergeoemtric differential equation itself, but rather by applying the
Beta method to Kummer’s formula .

Via (2.6):

$\frac{3^{1/2}2^{2/3}[\Gamma(1/3)]^{3}}{18\pi}=\frac{\Gamma(7/6)}{\Gamma(1/2)\cdot\Gamma(2/3)}\int_{0}^{1}\frac{x^{-1/3}(1-x)^{-1/2}}{1+x/8}dx$

We obtain the simplification:

$\frac{4\cdot
3^{1/2}\cdot\pi}{9}=\int_{0}^{1}\frac{x^{-1/3}(1-x)^{-1/2}}{1+x/8}dx$ (4.5)

Letting

$\displaystyle k=1,s=1,a=-1/3,b=-1/2,P(x)=(1+x/8)$ $\displaystyle
Q(x)=(8x-72),z=-72,a_{1}=8,a_{0}=-72$

We arrive at the identity:

$\int_{0}^{1}\frac{x^{-1/3}(1-x)^{-1/2}}{1+x/8}dx=\int_{0}^{1}\frac{x^{-1/3}(1-x)^{-1/2}(8x-72)}{-72-x(1-x)}dx$

And obtain (1.2).

##  5 Higher-Order Constants

There series we’ve been deriving so far involve constants that can be
expressed though the hypergeometric function of the form $2F1\\{1,x;y;z\\}$.
Such constants involve logarithms, root extractions, the gamma function, and
inverse trigonometric functions. Higher order constants such as Catalans’s
Constant, Zeta 3, etc can’t be expressed through $2F1\\{1,x;y;z\\}$, but
require a more generalized hypergeometric function of the form:
${}_{q+1}F_{q}\left\\{1,x_{1},\ldots,x_{q};y_{1};\ldots;y_{q};z\right\\}$.

Unfortunately, the algorithm (1.4) isn’t as efficient at generating rapidly
converging series for constants that are only expressible though a
hypergeometric function where $q>1$. To the best of my knowledge for (1.4) to
work on higher order constants such as Zeta 3, the value of ’ $s$ ’ must be
set to zero. Accelerating the convergence of a q>1 hypergeometric function via
(1.4) by setting $\mathrm{s}=0$ equates to taking multiple terms of

${}_{q+1}F_{q}\left\\{1,x_{1},\ldots,x_{q};y_{1};\ldots;y_{q};z\right\\}$ ’
$\mathrm{k}$ ’ terms at a time. So if $\mathrm{k}=2$, the series generated by
(1.4) would converge at a rate of $z^{2}$.

In this next example we will show how to accelerate the convergence of a
hypergeometric function of the form
$3F2\left\\{1,x_{1},x_{2};y_{1};y_{2};z\right\\}$ using (1.4), and then apply
it to derive faster converging infinite series for Catalan’s Constant and
$\Gamma(1/4)$.

Begin with the integrals:

$\displaystyle 3F2\left\\{1,x_{1},x_{2};y_{1};y_{2};z\right\\}=$ (5.1)
$\displaystyle\frac{\Gamma\left(y_{2}\right)}{\Gamma\left(x_{2}\right)\Gamma\left(y_{2}-x_{2}\right)}\int_{0}^{1}x^{x_{2}-1}(1-x)^{y_{2}-x_{2}-1}2F1\left\\{1,x_{1};y_{1};xz\right\\}dx$

$2F1\left\\{1,x_{1};y_{1};xz\right\\}=\frac{\Gamma\left(y_{1}\right)}{\Gamma\left(x_{1}\right)\Gamma\left(y_{1}-x_{1}\right)}\int_{0}^{1}\frac{y^{x_{1}-1}(1-y)^{y_{1}-x_{1}-1}}{1-yxz}dy$
(5.2)

Using the integral (5.2), choosing the values below for $\mathrm{k}$ and
$\mathrm{s}$, and performing polynomial long division we plug the following
into (1.4):

$\displaystyle
s=0,k=2,b=\left(y_{1}-x_{1}-1\right),a=\left(x_{1}-1\right),z=\frac{1}{x^{2}z^{2}}$
$\displaystyle
P(x,y,z)=(1-xyz),Q(x,y,z)=\left(\frac{y}{xz}+\frac{1}{x^{2}z^{2}}\right)$

And obtain:

$\displaystyle\int_{0}^{1}\frac{y^{x_{1}-1}(1-y)^{y_{1}-x_{1}-1}}{P(x,y,z)}dy=\int_{0}^{1}\frac{y^{x_{1}-1}(1-y)^{y_{1}-x_{1}-1}Q(x,y,z)}{\frac{1}{x^{2}z^{2}}-y^{2}}dy=$
(5.3)
$\displaystyle\frac{x^{2}z^{2}\cdot\Gamma\left(x_{1}\right)\Gamma\left(y_{1}-x_{1}\right)}{\Gamma\left(y_{1}\right)}\sum_{n=0}^{\infty}\frac{\left(x_{1}\right)_{2n}(xz)^{2n}}{\left(y_{1}\right)_{2n}}\left(\frac{x_{1}+2n}{xz\left(y_{1}+2n\right)}+\frac{1}{x^{2}z^{2}}\right)$

Plugging (5.3) into (5.2), and then plugging that result into (5.1) and
simplifying we get:

$\displaystyle\frac{\Gamma\left(y_{2}\right)}{\Gamma\left(x_{2}\right)\Gamma\left(y_{2}-x_{2}\right)}\int_{0}^{1}x^{x_{2}-1}(1-x)^{y_{2}-x_{2}-1}\sum_{n=0}^{\infty}\frac{\left(x_{1}\right)_{2n}(xz)^{2n}}{\left(y_{1}\right)_{2n}}\left(\frac{xz\left(x_{1}+2n\right)}{y_{1}+2n}+1\right)dx$
(5.4) $\displaystyle=3F2\left\\{1,x_{1},x_{2};y_{1};y_{2};z\right\\}$

Then we extract and convert the following integrals from (5.4) into pochhammer
& gamma notation via the beta function:

$\displaystyle\int_{0}^{1}x^{2n+x_{2}-1}(1-x)^{y_{2}-x_{2}-1}dx=\frac{\Gamma\left(2n+x_{2}\right)\Gamma\left(y_{2}-x_{2}\right)}{\Gamma\left(2n+y_{2}\right)}=\frac{\left(x_{2}\right)_{2n}\Gamma\left(x_{2}\right)\Gamma\left(y_{2}-x_{2}\right)}{\left(y_{2}\right)_{2n}\Gamma\left(y_{2}\right)}$
(5.5)

Then plugging (5.5) back into (5.4) and simplifying gives the final result:

$3F2\left\\{1,x_{1},x_{2};y_{1};y_{2};z\right\\}=\sum_{n=0}^{\infty}\frac{\left(x_{1}\right)_{2n}\left(x_{2}\right)_{2n}(z)^{2n}}{\left(y_{1}\right)_{2n}\left(y_{2}\right)_{2n}}\left(\frac{z\left(x_{1}+2n\right)\left(x_{2}+2n\right)}{\left(y_{1}+2n\right)\left(y_{2}+2n\right)}+1\right)$
(5.6)

The procedure used to derive (5.6) can be generalized:

${}_{q+1}F_{q}\left\\{1,x_{1},\ldots,x_{q};y_{1};\ldots;y_{q};z\right\\}=\sum_{n=1}^{\infty}\left(\left(\left(z\prod_{g=1}^{q}\frac{x_{g}+2n}{y_{g}+2n}\right)+1\right)\prod_{g=1}^{q}\frac{\left(x_{g}\right)_{2n}}{\left(y_{g}\right)_{2n}}\right)$
(5.7)

With regards to Catalan’s Constant it is well-known that

$3F2\left\\{1,1,\frac{1}{2};\frac{3}{2};\frac{3}{2};\frac{1}{4}\right\\}=\sum_{n=0}^{\infty}\frac{1}{(2n+1)^{2}\left(\begin{array}[]{l}2n\\\
n\end{array}\right)}=\frac{\pi}{3}\ln(2-\sqrt{3})+\frac{8}{3}\sum_{n=0}^{\infty}\frac{(-1)^{n}}{(2n+1)^{2}}$
(5.8)

Plugging the hyqpergeometric terms of (5.8) into (5.6) and simplifying gives
an infinite series for (5.8) that converges twice as fast:

$\sum_{n=0}^{\infty}\frac{1}{(2n+1)^{2}\left(\begin{array}[]{l}2n\\\
n\end{array}\right)}=\sum_{n=0}^{\infty}\frac{40n^{2}+54n+19}{((4n+1)(4n+3))^{2}\left(\begin{array}[]{l}4n\\\
2n\end{array}\right)}$ (5.9)

Using a transformation on (5.2) letting $\mathrm{k}=3$ gives:

$4\sum_{n=0}^{\infty}\frac{1}{(2n+1)^{2}\left(\begin{array}[]{l}2n\\\
n\end{array}\right)}=\sum_{n=0}^{\infty}\frac{6804n^{4}+17172n^{3}+15903n^{2}+6405n+956}{((6n+1)(6n+3)(6n+5))^{2}\left(\begin{array}[]{l}6n\\\
3n\end{array}\right)}$ (5.10)

In addition, formulas for $\Gamma(1/4)$ can be derived using (5.6) based on
accelerating a quadratic transformation. Begin with the following
quadratically transformed hypergeometric functions:

$2F1\left\\{\frac{1}{4},\frac{1}{4};1;\frac{-1}{8}\right\\}=\frac{\sqrt{\pi}}{2^{1/4}[\Gamma(3/4)]^{2}}$
(5.11)
$2F1\left\\{\frac{1}{4},\frac{3}{4};1;\frac{1}{9}\right\\}=\frac{\sqrt{3\pi}[\Gamma(1/4)]^{2}}{4\pi^{2}}$
(5.12)

(Note: Unlike equation (4.4), (5.11) and (5.12) can’t be derived with (1.4).
Proofs regarding quadratic hypergeometric transformations is given in [7].)

However, (5.11) is equal to
$3F2\left\\{1,\frac{1}{4},\frac{1}{4};1;1;\frac{-1}{8}\right\\}$ and (5.12)
can be written as
$3F2\left\\{1,\frac{1}{4},\frac{3}{4};1;1;\frac{1}{9}\right\\}$. Plugging
these values into (5.6) and simplifying via (2.5) gives

$\frac{36\sqrt{3}[\Gamma(1/4)]^{2}}{\pi^{3/2}}=\sum_{n=0}^{\infty}\frac{\left(640n^{2}+608n+147\right)(8n)!}{[(2n+1)!]^{2}(4n)!24^{4n}}$
(5.13)

$\frac{128\pi^{1/2}}{2^{1/4}[\Gamma(3/4)]^{2}}=\sum_{n=0}^{\infty}\left(448n^{2}+448n+127\right)\left(\frac{(1/8)_{n}(5/8)_{n}}{2^{n}(2n+1)!}\right)^{2}$
(5.14)

## References:

[1] Boris Gour’evitch, Jes’us Guillera Goyanes, Construction Of Binomial Sums
For pi And Polylogarithmic Constants Inspired By BBP Formulas (2006).

[2] D. Bailey, P. Borwein and S. Plouffe, On The Rapid Computation of Various
Polylogarithmic Constants, Math. Comp., 66(1997), 903-913.

[3] Travis Sherman, Summation of Glaisher- and Apéry-like Series, University
of Arizona (2000),http://math.arizona.edu/ ura/001/sherman.travis/series.pdf

[4] G. Almkvist, C. Krattenthaler and J. Petersson, Some new formulas for _,
preprint, Matematiska Institutionen, Lunds Universitet, Sweden. and Institut
fur Mathematik der Universitat Wien, Austria (2001),
http://www.arxiv.org/pdf/math.NT/0110238

[5] "List of hypergeometric identities" Wikipedia, The Free Encyclopedia.
Wikimedia Foundation, Inc. 10 Aug. 2004.

http://en.wikipedia.org/wiki/List_of_hypergeometric_identities

[6] Fabrice Bellard, _ page, http://fabrice.bellard.free.fr/pi/.

[7] Earl D. Rainville, Ph.D. "Special Functions", page 65, The MacMillan
Company, NY., 1960.

[8] Hakimoglu, Cetin, "An Algorithm for the Derivation of Rapidly Converging
Infinite Series for Universal Mathematical Constants" (2009,2021). Available
at SSRN: https://ssrn.com/abstract=3919892

[9] John M. Campbell et al., "On Some Series Involving the Binomial
Coefficients", 2019, https://arxiv.org/pdf/2306.16889.pdf

[10] Wenchang Chu, "DOUGALL’S BILATERAL 2H2-SERIES AND RAMANUJAN-LIKE
FORMULAE", 2011,
https://www.ams.org/journals/mcom/2011-80-276/S0025-5718-2011-02474-9/S0025-5718-2011-02474-9.pdf

[11] Wenchang Chu, "Ramanujan-like formulae for and via Gould–Hsu inverse
series relations", 2011,
https://link.springer.com/article/10.1007/s11139-020-00337-z

[12] Zhi-Wei Sun, "New series involving binomial coefficients",
2023, https://arxiv.org/abs/2307.03086

[13] K. A. DRIVER, "AN INTEGRAL REPRESENTATION OF SOME HYPERGEOMETRIC
FUNCTIONS", 2006,
https://www.emis.de/journals/ETNA/vol.25.2006/pp115-120.dir/pp115-120.pdf

[14] Necdet Batir, "On the series", 2005,
https://arxiv.org/pdf/math/0512310.pdf


# Wide Binary Orbits are Preferentially Aligned with the Orbits of Small
Planets, but Probably Not Hot Jupiters

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  Juliette Becker Department of Astronomy, University of Wisconsin-Madison, Madison, WI 53706, USA  Adam L. Kraus Department of Astronomy, University of Texas at Austin, Austin, Texas 78712, USA  Logan Pearce Steward Observatory, University of Arizona, Tucson, AZ 85721, USA  Karen A. Collins Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Malena Rice  Department of Astronomy, Yale University, New Haven, CT 06511, USA  Eric L. N. Jensen Department of Physics & Astronomy, Swarthmore College, Swarthmore PA 19081, USA  David Baker Physics Department, Austin College, Sherman, TX 75090, USA  Paul Benni  Acton Sky Portal (private observatory), Acton, MA USA  Allyson Bieryla Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Avraham Binnenfeld  Porter School of the Environment and Earth Sciences, Raymond and Beverly Sackler Faculty of Exact Sciences, Tel Aviv University, Tel Aviv, 6997801, Israel  Kevin I. Collins George Mason University, 4400 University Drive, Fairfax, VA, 22030 USA  Dennis M. Conti American Association of Variable Star Observers, 185 Alewife Brook Parkway, Suite 410, Cambridge, MA 02138, USA  Phil Evans  El Sauce Observatory, Coquimbo Province, Chile  Eric Girardin  Grand Pra Observatory, Switzerland  Joao Gregorio  Crow Observatory, Portalegre, Portugal  Tsevi Mazeh  School of Physics and Astronomy, Tel Aviv University, Tel Aviv, 6997801, Israel  Felipe Murgas  Instituto de Astrofísica de Canarias (IAC), E-38205 La Laguna, Tenerife, Spain  Departamento de Astrofísica, Universidad de La Laguna (ULL), E-38206 La Laguna, Tenerife, Spain  Aviad Panahi School of Physics and Astronomy, Tel Aviv University, Tel Aviv, 6997801, Israel  Francisco J. Pozuelos Space Sciences, Technologies and Astrophysics Research (STAR) Institute, Université de Liège, 19C Allée du 6 Août, 4000 Liège, Belgium  Astrobiology Research Unit, Université de Liège, 19C Allée du 6 Août, 4000 Liège, Belgium  Howard M. Relles  Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Fabian Rodriguez Frustaglia  American Association of Variable Star Observers, 49 Bay State Road, Cambridge, MA 02138, USA  Richard P. Schwarz Center for Astrophysics | Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Gregor Srdoc  Kotizarovci Observatory, Sarsoni 90, 51216 Viskovo, Croatia  Chris Stockdale Hazelwood Observatory, Australia  Thiam-Guan Tan Perth Exoplanet Survey Telescope, Perth, Western Australia  Curtin Institute of Radio Astronomy, Curtin University, Bentley, Western Australia 6102  William C. Waalkes Department of Physics and Astronomy, Dartmouth College, Hanover NH 03755, USA  Gavin Wang Department of Physics & Astronomy, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA  Justin Wittrock George Mason University, 4400 University Drive, Fairfax, VA, 22030 USA  Shay Zucker  Porter School of the Environment and Earth Sciences, Tel Aviv University, Tel Aviv 6997801, Israel

###### Abstract

Studying the relative orientations of the orbits of exoplanets and wide-
orbiting binary companions (semimajor axis greater than 100 AU) can shed light
on how planets form and evolve in binary systems. Previous observations by
multiple groups discovered a possible alignment between the orbits of visual
binaries and the exoplanets that reside in them. In this study, using data
from Gaia DR3 and TESS, we confirm the existence of an alignment between the
orbits of small planets $(R<6R_{\oplus})$ and binary systems with semimajor
axes below 700 AU ($p=10^{-6}$). However, we find no statistical evidence for
alignment between planet and binary orbits for binary semimajor axes greater
than 700 AU, and no evidence for alignment of large, closely-orbiting planets
(mostly hot Jupiters) and binaries at any separation. The lack of orbital
alignment between our large planet sample and their binary companions appears
significantly different from our small planet sample, even taking into account
selection effects. Therefore, we conclude that any alignment between wide-
binaries and our sample of large planets (predominantly hot Jupiters) is
probably not as strong as what we observe for small planets in binaries with
semimajor axes less than 700 AU. The difference in the alignment distribution
of hot Jupiters and smaller planets may be attributed to the unique
evolutionary mechanisms occuring in systems that form hot Jupiters, including
potentially destabilizing secular resonances that onset as the protoplanetary
disk dissipates and high-eccentricity migration occurring after the disk is
gone.

planet-star interactions

##  1 Introduction

Most stars in our galaxy reside in binary or higher-order stellar systems
(Raghavan et al., 2010), and many of these systems are known to host
exoplanets (Mugrauer, 2019). Using a variety of methods, astronomers have
detected exoplanets on a range of different orbits in binary systems, which
can generally be categorized into two main configurations: circumbinary, where
the planet orbits both stars in the binary system (e.g. Doyle et al., 2011;
Orosz et al., 2012), and circumstellar, where the planet orbits only one host
star (e.g. Anglada-Escudé et al., 2016; Mugrauer, 2019). In this paper, we
focus on exoplanets in circumstellar configurations in wide binary systems,
which happens to be the most common scenario for known exoplanets (Mugrauer,
2019). For the purposes of this paper, we define a wide binary to be any
binary star system with semimajor axis $\gtrapprox 100$ AU.

The effect of wide-binaries on the orbits of exoplanets in those systems has
been the subject of intense study for years. Generally, it is understood that
gravitational interactions with wide stellar binary companions can influence
the orbits of planets on long timescales. For example, the Lidov-Kozai-Von-
Zeipel (LKZ) mechanism (von Zeipel, 1910; Lidov, 1962; Kozai, 1962), induced
by the binary companion, could provide an explanation for at least some of the
orbits of hot Jupiters (Wu & Murray, 2003; Fabrycky & Tremaine, 2007;
Petrovich, 2015b; Dawson & Johnson, 2018), although other mechanisms that can
explain the existence of hot Jupiters exist, and it is not possible for the
LKZ effect to explain the existence of all Hot Jupiters (Becker et al., 2015;
Ngo et al., 2016). Moreover, occurrence rate studies from the Kepler mission
have shown that the presence of binary companions with separations less than
about 50 astronomical units can suppress planet formation (Kraus et al., 2016;
Moe & Kratter, 2021), presumably due to disruption of the protoplanetary disk.
And finally, theoretical work (Batygin, 2012; Lai, 2014; Zanazzi & Lai, 2018)
(and tentative observations, Jensen et al. 2020; Hjorth et al. 2021; Offner et
al. 2023) has suggested that binary companions can align protoplanetary disks
with the binary’s orbital plane, giving rise to systems with angular momenta
well aligned between planets and binary stars, yet misaligned with the stellar
spin axis.

Recently, several studies have focused on exploring this phenomenon in more
detail, and found evidence suggesting that the orbits of planetary systems and
the binaries they reside in are aligned111 It is important to distinguish
orbit-orbit alignment from spin-orbit alignment. Spin-orbit alignment is the
alignment between the stellar spin axis and exoplanet orbit, (e.g. Albrecht et
al., 2022). Orbit-orbit alignment is an alignment between the exoplanet orbit
and binary system orbit. (Behmard et al., 2022; Christian et al., 2022; Dupuy
et al., 2022). Christian et al. (2022) and Behmard et al. (2022) independently
found an alignment in TESS exoplanets and visual binaries identified from Gaia
DR2 data extending from 100 to over $10^{4}$ AU. Meanwhile Dupuy et al.
(2022), and more recently Lester et al. (2023), found a similar pattern
through astrometric monitoring of binary systems hosting planets discovered by
the Kepler mission. Thanks to the high resolution adaptive optics and speckle
imaging observations, these binaries could be resolved at closer separations
than Christian et al. (2022) and Behmard et al. (2022). These works all
concluded that closely orbiting planets and wide binaries show at least some
preferential alignment, and Behmard et al. (2022)’s results hinted at a
dichotomy between the orbit/orbit alignment states for small planets and hot
Jupiters.

These studies were limited in the conclusions they could draw, primarily
because any measurement of transiting planet orbit/binary orbit misalignment
must be statistical in nature. In particular, this is because for transiting
planet systems, it is only possible to measure the relative inclination
difference between binary and exoplanet orbit. Because inclinations are
measured relative to our line of sight, the difference in inclination between
the planet and stellar orbits is equal to the minimum misalignment between
orbits, so measuring an inclination difference of zero degrees, for instance,
could indicate a completely aligned or completely misaligned system. But by
measuring the inclination differences for a large sample of systems, a
preponderance of systems with small inclination differences implies a
preponderance of aligned systems. As a result of the ambiguity, large samples
are critical for confidently assessing the presence of aligned populations and
identifying any correlations with alignment or misalignment.

In this work, we therefore expand upon earlier works examining orbit/orbit
alignment with a significantly larger sample of planetary systems than were
previously studied. In particular, we perform a similar analysis to that of
Christian et al. (2022) and constrain the orbital inclination of a large
sample of wide binary star orbits in systems that have transiting planet
candidates from the TESS mission. Thanks to several years of new planet
candidate discoveries from TESS, new wide binaries detected in Gaia DR3, and
more precise astrometric solutions, we are able to draw stronger and more
nuanced conclusions on planet/binary alignment. While we find strong evidence
for alignment between small planets $(R<6R_{\oplus})$ and relatively close
binaries ($a<700$AU), we find no compelling evidence for alignment in systems
with giant planets or binary separations greater than $a\gtrsim 700$AU.

Our paper is organized as follows. We describe the construction of our samples
(both a sample of transiting planetary systems with wide binary companions,
and a control sample of wide binaries without any preference for transiting
planets) in Section 2 and describe our analysis of the data in Section 3. We
report our major results in Section 4 and discuss how they compare with
previous studies and what these trends imply for the early dynamical history
of these systems in Section 5. We then conclude in Section 6 with a summary of
our results and the most significant implications.

##  2 Observations/Data

###  2.1 Selection of wide binaries

We start by gathering the sample of wide binary stars identified using Gaia
EDR3 data by El-Badry et al. (2021). This updated sample contains
approximately one million common-proper-motion wide binary pairs within about
one kiloparsec of the Sun. This work builds upon a previous catalog compiled
by El-Badry & Rix (2018) (and used by Christian et al. 2022) using data from
Gaia DR2, which focused on systems within 200 pc, and included only about
50,000 binaries. Like Christian et al. (2022), we removed systems where the
astrometric solutions showed evidence of excess scatter, in particular by
removing systems where at least one star – either primary or companion – has a
Renormalized Unit Weight Error (RUWE) greater than 1.4, following the
recommendations of Lindegren (2018). We also removed putative binary systems
that might instead be chance alignments. We used the pseudo-probability of
chance alignments, $R$, as reported by El-Badry et al. (2021), to identify
these systems. We impose a relatively stringent requirement that the binary
systems have $R$ less than 0.01 to remain in our sample to ensure high purity.

We removed any systems with white dwarf companions. Post main sequence
evolution likely changes the architecture of the binary systems significantly
(Veras, 2016), and might affect the resulting inclination distribution. Using
only main sequence and red giant stars allows us to study potential primordial
or early-lifetime mechanisms that contributed to the observed alignment
distribution.

###  2.2 Selection of exoplanets

Like Christian et al. (2022), we source transiting exoplanet candidates from
the TESS Object of Interest catalog (Guerrero et al., 2021), which we
downloaded on 1 June 2022. The TESS exoplanets are listed by their identifiers
in the the TESS Input Catalog (TIC), which is based on Gaia DR2. We
crossmatched the Gaia DR2-based TESS planet candidates with the Gaia DR3-based
El-Badry et al. (2021) sample of wide binaries using the Gaia archive’s Gaia
DR2 to Gaia DR3 crossmatch.

After identifying the transiting planet candidates in binary systems via this
crossmatch, we excluded known and suspected false positives from our sample.
In particular, we performed the following steps:

  1. 1.

We queried the status of the planet candidates on ExoFOP-TESS to identify
which have been shown to be false positives or false alarms. These are
confidently known false positives, based on analysis of TESS data performed
after the planet candidates were initially identified, or follow-up
observations (mostly photometric or spectroscopic).

  2. 2.

Beyond the planet candidates that were identified as confident false positives
(FPs) on ExoFOP-TESS, we removed suspected FPs based on photometric follow-up
observations by Sub-group 1 (SG1) the TESS Follow-up Working Group (TFOP WG),
following Christian et al. (2022), Section 2.1.3.

The resulting sample contains 398 exoplanet candidates in binary systems (in
addition to the previously described cuts on binary RUWE and false positive
probability as described in section 2.1). This is a significantly larger
sample than the 67 objects used by Christian et al. (2022). Of these 398
systems, 59 are known exoplanet systems discovered before the TESS mission,
and 33 were confirmed to be exoplanets by the TESS mission. As of April 2024,
150 are designated by SG1 as verified planet candidates, verified candidates
with possible contaminated aperture, or cleared planet candidates, indicating
that most potential sources of a false-positive detection have been ruled out,
but not to a degree to designate the planet as confirmed.

Overall, only 79 systems in our sample that were not previously known to be
confirmed planets lack significant follow-up observations from TFOP SG1.

Figure 1: Comparisons of planet radius and parallax distributions from our
sample to that of Christian et al. (2022). Left: a comparison of the
histograms of planet radius in the new and old samples. Right: a comparison of
the joint distribution of the parallax and planet size between the sample in
this work and that of Christian et al. (2022). The significant increase in new
giant planets in our sample is primarily due to the El-Badry et al. (2021)
catalog extending to further distances (smaller parallaxes) than that of El-
Badry & Rix (2018). Large transiting planets are inherently rarer, but easier
to detect around more distant/faint stars. In the right-hand panel, the points
corresponding to the same systems in the old and new samples do not precisely
line up because of the difference of Gaia DR2 and DR3 parallaxes and changes
in the TESS planetary radii as new data have been collected.

The increase in number of systems in our sample is accompanied by changes in
the characteristics of the typical system in our sample compared to that of
Christian et al. (2022). Figure 1 shows histograms that illustrate these
changes. Compared to the sample of Christian et al. (2022), ours is a) more
dominated by giant planets, especially hot Jupiters, and b) composed of more
distant stars. These changes are explained by the fact that the El-Badry et
al. (2021) catalog extends to greater distances than previous work, thereby
increasing the number of faint binary systems. Moreover, the significant
increases in the TOI catalog are largely due to the TESS faint star search
(Kunimoto et al., 2022), which increases the magnitude limit to which planet
candidates are reported from the MIT Quick Look Pipeline from 10.5 to 13.5 in
TESS-band magnitudes. Both of these changes increased the number of faint,
distant stars in our sample, and at these faint magnitudes, TESS is primarily
sensitive to short-period gas giant exoplanets. As a result, our sample has a
much greater proportion of hot Jupiter exoplanets than that of Christian et
al. (2022).

###  2.3 Control Sample

In order to control for unknown biases in the inclination distribution of
binaries from Gaia DR3, we build a control sample of systems without detected
exoplanets. As in Christian et al. (2022), we selected the members of our
control sample by searching for “similar” star systems to our exoplanet hosts.
We quantified the similarity between our planet-hosting systems and putative
members of our control sample using a slightly different metric $\mathcal{M}$
than Christian et al. (2022). In particular, we defined:

$\mathcal{M}=(\Delta G_{1})^{2}+(\Delta
G_{2})^{2}+\left(\frac{\Delta\varpi_{1}}{\varpi_{1,\rm
exo}}\right)^{2}+\left(\frac{\Delta s}{s_{\rm exo}}\right)^{2}$ (1)

where $G_{1}$ and $G_{2}$ are the Gaia G-band magnitudes of the primary and
secondary stars, $\varpi_{1}$ is parallax of the primary stars in each binary
system, and $s$ is binary angular separation. Here, we use $\Delta$ to refer
to the difference in each quantity between that of the exoplanet-host system
and the values of putative control-sample members, and the subscript exo
denotes the particular value of the exoplanet-hosting system. Matching on
apparent magnitude and parallax is a proxy for matching the luminosity of the
binary system. There is tentative evidence that the inclination distribution
of binary systems changes as a function of mass ratio (Gerbig et al., 2024),
so this is an important quantity to control for. The primary difference
between this metric and that used by Christian et al. (2022) is that we do not
compute the normalized difference in magnitude, as the magnitude already is
defined on a log-scale.

For each binary system in our sample of binaries with exoplanets, we select
the ten systems from the El-Badry et al. (2021) sample that, after passing the
initial quality cuts of RUWE and R described in Section 2.1, have the lowest
value of $\mathcal{M}$.

##  3 Analysis

###  3.1 Stellar Mass Determination

Figure 2: Histogram of the primary companion and secondary companion stellar
masses. The primary companion of a binary star system is defined as the star
which the exoplanet orbits.

First, we estimated the masses of each stellar component (primary and
secondary) in both our exoplanet and control samples. Like Christian et al.
(2022), we use two different methods: an empirical relation for low-mass
M-dwarfs, and fitting to models for higher-mass stars.

For the low-mass M-dwarf stars – that is, stars whose absolute $K$-band
magnitudes are between 4.5 and 10.5 – we followed Christian et al. (2022) and
estimated masses using the empirical $M_{\odot}-K$ relationship of Mann et al.
(2019), implemented in an accompanying python
package222https://github.com/awmann/M_-M_K-. For the absolute magnitude
determination for M-dwarfs, we use the geometric distances of Bailer-Jones et
al. (2021).

For the higher-mass stars, we again followed Christian et al. (2022) and
estimated masses using the `Isochrones` package (Morton, 2015), which infers
fundamental stellar parameters using nested sampling to fit stellar
evolutionary and atmospheric models given data like broadband apparent
magnitudes, spectroscopic metallicity, and parallax. For each star in each of
our binary systems, we input the G,BP, and RP Gaia magnitudes, as well as J,H,
and K magnitudes from the 2MASS survey (if available), and the star’s Gaia
parallax. Unlike Christian et al. (2022), we include only metallacity
estimates from Gaia, and do not make use of any other spectra (either from
ground-based surveys or follow-up observations). Gaia DR3 provides low-
resolution spectroscopy for certain systems, including derived metallicity.
When available, we use the provided spectroscopic metallicity values from Gaia
DR3 and cut off values below -1 (Fouesneau et al., 2022). We find for the
stars in our sample with both spectroscopic and photometric metallicities, the
values generally agree, and choose to use the spectroscopic metallicities.

The mass distribution of the primary and secondary stars is plotted in Figure
2. A large fraction of the secondary stars are M-dwarfs, while the primary
stars generally have masses above one stellar mass.

###  3.2 LOFTI modeling

Like Christian et al. (2022), once we had mass estimates for all of the stars
in both our exoplanet and control samples, we we fitted the orbits of the
binary systems using Linear Orbits for The Impatient (Pearce et al., 2020,
LOFTI). LOFTI uses a Monte Carlo approach to perform rejection sampling in a
simplified 4D orbital parameter space (Blunt et al., 2017). For each system,
given the stellar masses, LOFTI randomly generates model orbits, compares them
to observational data (i.e. positions, parallax, proper motions, and radial
velocities), and either accepts or rejects the samples based on the the
goodness of fit. LOFTI repeats this process until a user-specified number of
accepted samples is reached; the distribution of the accepted samples then
represents the posterior probability distribution for orbital parameters.

We ran LOFTI on each system in both our exoplanet and control samples. For
each system, we input the positions, proper motions, parallax, and radial
velocities (where applicable) for each star from Gaia DR3. Because of the
significantly larger number of systems in our sample compared to Christian et
al. (2022), we chose to decrease the number of accepted samples per system
from 100,000 to 10,000. We found that 10,000 orbits were plenty for us to make
reasonable statistical estimates, and this sped up computation by a factor of
10.

For a small number of systems (3 in our exoplanet sample, and 61 in our
control sample), LOFTI was unable to generate a sufficient number of accepted
samples (typically less than 5 accepted samples per 24 hours of runtime on a
single CPU core). Often, when LOFTI is unable to find good orbital solutions
like it, it is a sign that the system is not a real, physically bound binary,
and may instead just be a chance alignment. Given the small number of these
systems, the computational expense that would required to reach 10,000 orbits,
and the fact that the systems are likely not physically bound anyway, we chose
to exclude them from our final samples.

Figure 3: Histogram of the minimum inclination difference between binary
orbits and the line of sight ($|90^{\circ}-i|$) for both a sample of binaries
with exoplanets and a control sample without known exoplanets. An istropic
distribution of inclinations generated using $p(i)\sim\cos(i)$ is also plotted
in red. These histograms show our full sample, regardless of the planet radius
and binary semimajor axis. Unlike Christian et al. (2022), we see no clear
evidence for preferential alignment in our full sample. We attribute this to
the fact that the Christian et al. (2022) sample was primarily composed of
small planets, which do show alignment (see Section 4.2.1) while our full
sample has many more hot Jupiters, which do not show a preferential alignment.

##  4 Results

Figure 4: Four histograms of inclination binned by semimajor axis (left versus
right) and radius (top versus bottom). The x-axis of each histogram is
$|90-i|$. Orange represents the control sample, while blue represents the
sample of binaries with exoplanets. The background color of each plot
represents the p-value from the K-S test performed between control and the
sample with exoplanets. We see strong statistical evidence for alignment in
small planets in binary systems with small semimajor axes, but in none of the
other regions.

###  4.1 No Clear Preferential Alignment between Binary and Planet Orbits in
the Full Sample

Once we had completed our orbital modeling for each system in both the
exoplanet sample and the control sample, we calculated $|90^{\circ}-i|$, the
minimum misalignment between the orbit of a transiting exoplanet in each
system and the binary orbit. In general, if there is alignment between
planetary orbits and binary orbits, we would expect to find an overabundance
of systems with $|90^{\circ}-i|$ near $0^{\circ}$ degrees for the exoplanet
sample. We compared the distributions of $|90^{\circ}-i|$ for the exoplanet
sample and the control sample, and show the resulting histograms in Figure 3.
We also plot an isotropic distributions of inclinations.

While there is a hint of an overabundance of exoplanet-hosting systems near
$|90^{\circ}-i|=0^{\circ}$, the effect is considerably weaker than the
preferential alignment seen by Christian et al. (2022) and is not
statistically significant. A Kolmogorov-Smirnov test is unable to confidently
reject the null hypothesis that the exoplanet sample’s inclinations are drawn
from a different population than the control sample’s inclination $(p=0.067)$.
Both the control sample and the sample with exoplanets are significantly
different from an isotropic distribution. This suggests that there are some
underlying selection biases in our sample of binaries, so any statistical
statement on inclination alignment should be made with respect to the control
sample, not an isotropic distribution.

###  4.2 Analysis of Sample Split in Planet Radius and Binary Semimajor Axis

This lack of an overabundance of exoplanets with $|90^{\circ}-i|\approx
0^{\circ}$ is somewhat surprising given the evidence from multiple independent
datasets (Christian et al., 2022; Dupuy et al., 2022; Behmard et al., 2022;
Lester et al., 2023) of preferential alignment between planet and binary
orbits in at least some parts of parameter space. We therefore investigated
whether subsets of our sample do show a preferential alignment. In particular,
Christian et al. (2022) noticed that the alignment they observed was strongest
for closely orbiting binaries (with semimajor axes less than about 700 AU),
while Behmard et al. (2022) found tentative evidence for misalignment of gas
giant planets.

We therefore chose to divide our sample along these two axes: planet radius
(greater than or less than 6 $R_{\oplus}$) and binary semimajor axis (greater
than or less than 700 AU). The divide of 6 $R_{\oplus}$ is chosen to roughly
divide small and giant planets. We explore the sensitivity to the cutoff of
semimajor axis cutoff in section 4.3.1.

The result of this is four subsamples, each containing between 30 and 86
objects. In a small number of cases, multi-planet systems host at least one
planet smaller than 6 $R_{\oplus}$ and at least one planet larger than 6
$R_{\oplus}$. In these cases, we chose to assign the system to both sub-
samples. While this may raise concerns about double-counting, each of the sub-
samples is analyzed independently, and we also confirmed that this strategy
did not significantly affect any of the results in our paper.

Figure 4 shows histograms (similar to that shown in Figure 3) comparing the
minimum misalignment between binary orbits and the line of sight for the
control sample and for our exoplanet sample in each of the four sub-samples.

We calculate the Kolmogorov-Smirnov (KS) statistic for each of the four sub-
samples and visually represent the p-value in the plots by shading the
background color darker for lower p-values. No correction (e.g. Bonferroni
correction or any of its variants) was applied for doing multiple comparisons.
The derived p-values are $a<700$ AU, $r<6R_{\oplus}:4\times 10^{-6}$, $a<700$
AU, $r>6R_{\oplus}:0.4$, $a>700$ AU, $r<6R_{\oplus}$: $0.09$, $a>700$ AU,
$r>6R_{\oplus}$: 0.29.

The Anderson Darling test, a variant of the KS test which works better in some
cases where the difference is at the edges of distributions, returns p-values
of $a<700$ AU, $r<6R_{\oplus}:0.001$, $a<700$ AU, $r>6R_{\oplus}:0.2$, $a>700$
AU, $r<6R_{\oplus}$: $0.1$, $a>700$ AU, $r>6R_{\oplus}$: 0.2. For the
remainder of this paper, we will use the KS statistical test, since the
results are qualitatively similar when using the Anderson Darling test.

Figure 5: (a) Plot of median inclination vs. median semimajor axis for
binaries with exoplanets from Gaia DR3. Some exoplanets do not have reported
radii values so are represented as maroon squares. Since the errorbars on
inclination are the 16th and 84th percentiles of $|90-i|$, while the median
values are reported as $|90-\overline{i}|$ where $\overline{i}$ is the median
value of inclination, occasionally an errorbar will go under zero degrees in
which case we cut it off. (b) The same as (a), but for the control sample.
(c), (d) are KDEs of systems with planets that have with radii greater than or
less than $6R_{\oplus}$, respectively. (e) is a KDE of the control sample
systems in panel (b). The KDEs are drawn from 100 random samples from each
binary to reflect the underlying error present in the sampling procedure. The
precise difference in the distributions of various distributions is best seen
in Figure 6.

####  4.2.1 Strong Alignment in Small Planet Systems with Binary Semimajor
Axes Below 700 AU

Figure 4 shows a strong alignment between the orbits of small planets and
binary companions with semimajor axes less than about 700 AU (lower left
panel). The Kolmogorov-Smirnov test rejects the null hypothesis that the
inclinations of the control sample binaries are drawn from the same population
as the inclinations of close binary systems containing small planets (with
p-value of $1.2\times 10^{-6}$). This finding confirms the previous hints from
Christian et al. (2022) and Behmard et al. (2022), and increases the level of
statistical significance by several orders of magnitude.

The analysis shown in Figure 4 shows with high confidence that there is
preferential alignment between the orbits of small planets with semimajor axes
less than 700 AU and the orbits of the binary systems they reside in, but this
does not take into account the uncertainties in both the binary inclination
angle and semimajor axis, which in many cases are fairly large. We therefore
performed analysis to visualize the impact of uncertainties on these
parameters. We randomly sampled 100 points from the posterior probability
distributions for the binary inclination angles and semimajor axes for each
system in both the exoplanet sample and the control sample. We calculated a 2D
Kernel Density Estimate (KDE) of the mutual distribution of these parameters
for both samples, which we show in Figure 5. Subtracting the KDEs of the two
samples (Figure 5, panel e) reveals a significant excess in probability mass
at near edge-on binary orbits and semimajor axis less than 1000 AU. Evidently,
the pattern we noticed based on the distribution of median values of a given
parameter for each system still exists when taking uncertainties on these
parameters into account. We also see evidence of a dichotomy between small
planets and large planets when viewing the sample in this manner (see Section
4.3.1 for more details). Figure 6 shows the same KDEs, but now split between
small planets (see Figure 6 panel a, which show a strong excess in probability
mass for edge-on orbits at small semimajor axes) and large planets (see Figure
6, panel b, which shows no clear excess at edge-on orbits, but some preference
for misaligned orbits with respect to the control samples at semimajor axes
between 1000 and 10,000 AU). Finally, Figure 6, panel c shows the difference
between panels a and b, highlighting the overabundance of well-aligned systems
with small planets.

Figure 6: A similar arrangement to Figure 5 (c-e) except comparing differences
between populations. The relative density difference is calculated as the
difference in the normalized density of Figure 5 between two different
populations. (a) shows binaries with $r<6R_{\oplus}$ compared to the control
sample, (b) the same but with $r>6R_{\oplus}$, and (c) comparing
$r<6R_{\oplus}$ and $r>6R_{\oplus}$. For example, darker red in (a) indicates
relatively more small planets compared to the control sample, while lighter
colors mean that the small planet and control sample have similar densities.
The relative biases in the semimajor axis distributions of large and small
planets are apparent in panel (c). Because larger planets are found at further
distances, the semimajor axes of the binary systems they are found in is
systematically larger. Figure 7: The variation of the p-value between binaries
with exoplanets above and below a specific semimajor axis and the control
sample, as measured through a KS test. Explicitly, denoting the control sample
as $c$ and the binaries with exoplanets as $b$, the blue line is computed as
$\text{KS}(b[a_{i}<a\text{ and }r<6],c)$, the orange line as
$\text{KS}(b[a_{i}>a\text{ and }r<6],c)$, the red line as
$\text{KS}(b[a_{i}<a\text{ and }r>6],c)$, and the green line as
$\text{KS}(b[a_{i}>a\text{ and }r>6],c)$. Because we are doing many
comparisons, per Schlaufman (2015), the minimum of the p-value distribution
should not be interpreted as an exact probability, but the dip at $\approx$
700 AU suggests that there is some change in the inclination distribution
around roughly that semimajor axis.

####  4.2.2 No Clear Evidence for Alignment in the Largest-Separation Binary
Systems

Christian et al. (2022) found that the alignment they observed in their sample
was largely driven by systems with binary semimajor axes less than 700 AU. We
find a similar result. Neither of the sub-samples with binary semimajor axes
greater than 700 AU show statistically significant evidence for alignment. The
Kolmogorov-Smirnov test yields $p=0.29$ for the small planets and $p=0.09$ for
the large planets. Christian et al. (2022) interpreted the fact that alignment
is strongest at semimajor axes less than 700 AU as tentative evidence that
dynamical realignment of protoplanetary disks by the binary companions may
explain the pattern they observed. We reassess this hypothesis in light of our
increased sample size in Section 5.2.

We note that our nominal value of 700 AU for the cutoff binary semimajor axis
below which alignment becomes strong is not a precisely measured value.
Christian et al. (2022) identified this cutoff by visually inspecting a plot
of minimum alignment versus binary semimajor axis, and we use this value to
split our sample for the sake of consistency with their analysis. We wish to
emphasize that our choice to split the samples at 700 AU should not be taken
as evidence that there is a physical cutoff at precisely 700 AU, or that
indeed there are two distinct populations of systems separated by some
semimajor axis at all. Despite the increased sample size over our previous
paper, there are still only 40 small planets with semimajor axes below 700 AU,
so determining the exact location of the cutoff is difficult.

We performed an analysis inspired by Buchhave et al. (2014) to assess the most
likely location of a transition at which binaries are preferentially aligned
with short-period planets. We took all four of our subsamples and repeated our
analysis calculating the p-value of a Kolmogorov-Smirnov test between our
exoplanet and control sample when we split the samples at different binary
semimajor axis values. The p-values we measured are plotted against the
semimajor axis above which we excluded systems from our calculation in Figure
7. We indeed see a minimum in the p-value for small planets with small binary
separations around 700 AU, but the minimum is broad; from this, we infer that
the cutoff semimajor axis above which alignment disappears, if it exists, is
unlikely to be above $10^{3}$ AU. We note that the decrease in p-value for
small planets past $10^{2.5}$ AU is mainly due to an increasing sample size
rather than different sample characteristics at those semimajor axis values.
Per Schlaufman (2015), we caution against reading too much into the exact
minimum p-value in Figure 7, and we note that it is also entirely possible
that there is instead a gradient of systems from aligned to not aligned that
exists somewhere around 700 AU.

There is also a dip in the p-value when separating large planets above and
below 5000 AU. The p-value at this dip is less than 0.01, but given the many
comparisons taking place across the range of semimajor axis cutoffs, this
cannot be interpreted to be statistically significant.

###  4.3 No Clear Evidence for Anti-Alignment in Large Planet Systems

We inspected the histograms for large planets (with radii greater than 6
$R_{\oplus}$). Visually, there is an apparent overdensity in the number of
systems with semimajor axes greater than 700 AU with minimum inclination
differences of about 35 to 45 degrees. This feature is intriguing at first
glance because this is close to the maximum of the inclination range of Lidov-
Kozai oscillations (see e.g. Naoz 2016), but there is no evidence this is a
statistically significant feature. Indeed, when we instead plot the histogram
of $\sin(|90-i|)$, where the isotropic distribution takes on the form of a
uniform distribution, the excess does not not stand out above the noise in the
histogram, and no longer appears visually statistically significant.

####  4.3.1 No Clear Evidence for Alignment in Large Planet Systems

Finally, we see no clear evidence for alignment between wide binaries with
semimajor axes below 700 AU, and large planets in Figure 4. The Kolmogorov-
Smirnov test yields a p-value of 0.4 for this subsample. This is in contrast
to the strong alignment we see between the orbits of small planets and wide
binaries with semimajor axes less than 700 AU. The lack of a statistically
significant alignment for large planets at close semimajor axes could be the
result of either an astrophysical difference between small and large planet
systems, or it could be an artifact of the samples involved (either the large
planet sample is too small for adequate comparisons, or it is contaminated).
We therefore strove to determine which is the case.

First, we considered whether the samples are simply too small to make
definitive comparisons. The sample of large planets at small semimajor axes
contains only $\approx 30$ systems, and the sample of small planets at small
semimajor axes is not much larger ($\approx 40$). This makes it difficult to
directly compare the two distributions with high statistical confidence. We
tested whether the sample size is sufficient to confidently distinguish the
two inclination distributions by performing a KS test directly comparing the
large and small planet samples. On face value, the difference in distributions
of close-in large and small planets is border-line statistically significant,
with a KS test p-value of 0.012, implying that the two distributions do seem
significantly different, and that the samples are indeed large enough to draw
conclusions about differing degrees of alignment.

Second, we considered the more vexing possibility that the difference between
the two samples is caused by a different false positive rate between the two
samples, and therefore different levels of false positive contamination. The
vast majority of the large planets in our sample are hot or warm Jupiters,
which are intrinsically rare, and are difficult/impossible to distinguish from
other astrophysical scenarios (like eclipsing brown dwarfs and low-mass stars)
from TESS data alone. Therefore, the TESS false positive rate for these planet
candidates ($\approx$ 35% for hot Jupiters in particular) is considerably
higher than for small planet candidates ($\approx$ 15%; Zhou et al., 2019).
One mitigating factor for our study is that these rates are probably
overestimates, since we have already discarded a number of false positives
from our sample thanks to reconnaissance ground-based observations. 72 out of
250 of the hot Jupiters in our sample are confirmed planets, and more are
suggestively planets from follow-up observations. Since many planets are
already confirmed (or confirmed as false-positives), a 35% false-positive rate
is a high upper bound for the true false positive rate.

In an attempt to quantify the effects of the larger false positive rate of
large planets, we devise a bootstrap test based on the sample of small
planets. Our analysis seeks to answer the question: “If the inclination
distribution for large planets were drawn from the same inclination
distribution for small planets, would false positive contamination cause us to
see as large of a difference between the populations as we observe?” To do
this, we performed a test where we drew sample inclinations from the
distribution of inclinations we found for small planets, added in some
fraction of randomly chosen inclinations from our control sample, and
calculated the KS statistic between the new, simulated sample of giant planets
(with false positives) and the control sample. In recognition of the strong
dependence on the assumed false positive rate and the inexact nature of the
cutoff at about 700 AU in binary separation, we perform these simulations for
a variety of different values of these parameters.

In particular, we took the following steps:

  1. 1.

We selected the false positive rate ($X$) and the semimajor axis cutoff ($Y$)
for a particular bootstrap test. The procedure is computationally taxing, so
we computed the bootstrap procedure on a grid of only 20 false positive points
and 500 semimajor axis points.

  2. 2.

Given the chosen semimajor axis cutoff $Y$, we counted the number of systems
with planets larger than 6 $R_{\oplus}$ and semimajor axis less than $Y$. We
call this number of large-planet systems $A$.

  3. 3.

We then identify the sample of small planet systems with semimajor axes less
than $Y$.

  4. 4.

We randomly select $\lfloor A(1-X)\rfloor$ systems from the sample of small
planet identified in the previous step and record the inclinations of these
systems. Here, $\lfloor\rfloor$ denotes the “floor” operator since only an
integer number of samples can be used in a given comparison. These
predominantly well-aligned binary inclinations simulate a well-aligned
distribution of large planet systems.

  5. 5.

We then randomly select $\lceil AX\rceil$ systems from the control sample and
record the inclinations of these systems, adding them to the sample of
inclinations identified in the previous step, and bringing the total number of
inclinations in the synthetic sample to $A$. These randomly aligned systems
from the control sample simulate the effects of a population of false
positives with no preferential alignment among the hot Jupiters. Here,
$\lceil\rceil$ denotes the “ceiling” operator.

  6. 6.

We calculate the KS statistic comparing the simulated population of giant
planets (with false positives mixed in) and the control sample for giant
planets, and record the p-value. We repeat this test 10,000 times for each
value of $X$ and $Y$, and record the median p-value and fraction of the tests
where the p-value is below 0.05 for each $(X,Y)$ pair.

Figure 8 shows the results of these simulations for a variety of different
false positive rates and semimajor axis cutoffs. The left-hand panel shows the
median p-value for the suite of tests, and the right-hand panel shows the
fraction of the tests where the p-value is below 0.05. As expected, for low
false positive rates (where the fewest misaligned false positive systems are
mixed in) and small values of semimajor axis cutoff (where the alignment is
strongest), we see low median p-values and high fractions of the time yielding
a significant difference between the simulated sample and the control sample.
In particular, for the nominal values for these parameters (roughly 25% false
positive fraction and 700 AU cutoff), we expect that we should usually see a
statistically significant difference between the exoplanet sample and the
control sample. For injections with a cutoff of 700 AU and false positive
rates less than or equal to 25%, the minimum of percentage of samples in a
given injection with p-value less than 0.05 is 90%, a p-value of 0.1 is 96%, a
p-value of 0.2 is 98%, and a p-value of 0.3 is 99%. In reality, we see a KS
test p-value when we compare the exoplanet and control sample considerably
higher than these typical values (0.4), implying that false positive
contamination most likely cannot explain the difference between the small
planets and large planets. This suggests that large planets in binary systems
are probably not aligned to the same degree as small planets.

A 25% false-positive rate does not include the built-in false-positive rate
for small planets. If, say, the small planet rate is 15%, then this would
correspond to a 40% false-positive rate for hot Jupiters. As a reminder,
follow-up observations have already confirmed that some of the hot Jupiters in
our sample are planets, so a simulated false-positive rate of 25% is a
generous upper bound.

Figure 8: Results of simulations to determine whether higher false positive
rates for giant planets can explain the difference in the inclination
distribution compared to small planet systems (described in Section 4.3.1).
Left: an image showing the median p-value from 1000 simulations performed for
a given false positive fraction and semimajor axis cutoff. Right: an image
showing the percentage of p-values from each of the 1000 simulations that are
below 0.05. The horizontal and vertical lines visible in the figure are the
separations between the grid points used in the simulation. In both plots, a
blue x marks our expectation for the semimajor axis delineation and false-
positive rate. Our grid points have a much finer resolution in semimajor axis
than false positive rate. We find that for plausible values of the false
positive rate (expected to be roughly less than 25%) and semimajor axis cutoff
($\approx$ 700 AU), we would probably still expect to see a significant
alignment if giant planets were similar to small planets. Since we do not see
a significant alignment, we infer that giant planets are probably less
strongly aligned than small planets.

###  4.4 No dependence on temperature of the host star

Figure 9: Left: The distribution of $|90-i|$ for all systems with exoplanets
when the temperature of the host star of the exoplanet has a temperature
greater than or less than 6100 K, respectively. Unlike the distribution of
spin/orbit angles for hot Jupiters, we find no evidence of a difference in the
distribution of planet/binary orbit angles above and below the Kraft Break.
Right: The same as the left but restricted to less than 700 AU, and with
counts shown instead of relative density.

Finally, we tested whether the degree of alignment between planet and binary
orbits depends on the temperature of the host star in question. We are
motivated to perform this test because temperature of the host star is an
important predictor of whether the spin/orbit angle (obliquity) of hot Jupiter
systems is aligned or misaligned. Famously, hot stars above the Kraft break
($\sim 6100$ K) host hot Jupiters that span a wider range of obliquities than
cool stars (e.g. Winn et al., 2010; Schlaufman, 2010; Albrecht et al., 2012).
More recently, Rice et al. (2024) suggested that the orbits of close-in binary
star and exoplanet orbits should be fully aligned at small semimajor axes for
temperatures below the Kraft break, and less alignment should be seen above
the Kraft break. However, we find no dependence of the alignment between
planet orbits and wide binary orbits on the host star’s temperature, either
for our full sample, or only systems with large planets (dominated by hot
Jupiters). We plot the distribution of inclinations for host stars with
temperatures less than and greater than 6100 K in Figure 9. For all stars, a
KS test between host stars with temperatures less than and greater than 6100 K
returns a p-value of 0.55. For stars with large planets, the KS test gives a
p-value of 0.80 and for small planets, a p-value of 0.54. We attribute the
lack of a statistically significant difference here to the very small number
of hot stars with close-in binary companions ($a<700$ AU) in our sample. There
are only 9 systems in our sample with $T>6100$ K and $a<700$ AU.

##  5 Discussion

###  5.1 Comparison to Previous Work

The results of this work are largely in agreement with those of previous
papers (Christian et al., 2022; Dupuy et al., 2022; Behmard et al., 2022;
Lester et al., 2023), but strengthen and add some nuance to their conclusions.

Figure 10: Left:  Comparison of the distributions of binary inclination for
systems with small planets with binary semimajor axes less than 700 AU (Left)
and with binary semimajor axes greater than 700 AU (Right) between our sample
and that of Christian et al. (2022). In both panels, our results are
consistent with those of Christian et al. (2022), but our considerably larger
sample lets us make inferences with higher statistical confidence.

Previously, Christian et al. (2022), found that wide binaries (binary systems
with semimajor axes greater than 100 AU) and planets of all sizes showed
statistically significant alignment, and that the alignment was most prevalent
in systems with binary semimajor axes less than 700 AU. Now, thanks to a
considerably larger sample, with a greater diversity of planetary types, we
see the alignment identified by these previous works is predominantly in
systems containing small planets and binary systems with small semimajor axes.
This is in good agreement with the results of both Christian et al. (2022) and
Dupuy et al. (2022), who studied closer binary systems than Christian et al.
(2022) and found a higher degree of alignment than in the overall Christian et
al. (2022) sample. In Figure 10, we plot distributions of inclination for our
sample compared to Christian et al. (2022). We quantified the consistency
between our sample and that of Christian et al. (2022) by performing a KS
tests comparing the two samples. A KS test comparing the distribution of small
planets with $a<700$ AU in our sample versus that of Christian et al. (2022)
returns a p-value of 0.50, and a KS test comparing the samples for small
planets with $a>700$ AU yields a p-value of 0.59, indicating good consistency.

We also are able to corroborate the tentative evidence found by Behmard et al.
(2022) that hot Jupiter systems do not show the same alignment as small planet
systems, and in fact may be misaligned. With a larger sample than Behmard et
al. (2022), we also see no evidence for alignment for systems with large
planets, and see tentative evidence for preferential misalignment in Figure 6.
These findings should continue to strengthen as the planets discovered by TESS
in this sample continue to receive follow-up observations and mitigate the
risk that false positive contamination is affecting our results.

###  5.2 Implications for dynamical conclusions

Christian et al. (2022) suggested that a disk torquing mechanism, influenced
by energy dissipation in the protoplanetary disk due to a binary companion,
would result in alignment of the orbits of planets forming the disk with the
binary orbit. This seemed to explain the observed alignment cutoff for planets
of all masses at approximately 700 AU.

However, our new expanded data set reveals a contrast not apparent in our
previous smaller sample: while systems containing smaller planets exhibit
alignment within this range, systems containing large planets (predominantly
hot Jupiters333We note that while our large-planet sample nominally contains
both hot Jupiters and warm Jupiters, we treat it as dominated by hot Jupiters.
We attempted to test whether hot Jupiters and warm Jupiters showed different
behavior by separating them using the definition of Huang et al. (2016), but
there too few (only 23) warm Jupiter systems to draw definitive conclusions.
Adding warm Jupiter systems from Kepler and from ongoing TESS discoveries as
the survey baseline increases may enable such analysis in the future.) do not.
This suggests that systems containing small planets and those containing hot
Jupiters undergo different evolutionary paths. Two potential explanations
emerge: dynamical excitation of star-planet misalignments during the late
phases of the disk lifetime and disk dissipation, or subsequent dynamical
evolution of the hot Jupiter orbit over much longer timescales after the disk
has dissipated.

First, the picture described in Christian et al. (2022) of disk-binary orbits
being aligned due to energy dissipation does not account for the potentially
destabilizing effect of secular resonances, which may occur when the
fundamental frequencies mediating the interactions in the planetary system
evolve due to dissipating disks (Spalding & Batygin, 2014; Epstein-Martin et
al., 2022). When two precession frequencies in the system become commensurate
due to this evolution, it can lead to large-scale excitation of orbital
parameters. In particular, systems containing hot Jupiters, where a large
amount of angular momentum is concentrated at one orbital radius, may be
strongly affected by secular resonances as their parent disks dissipate
(Zanazzi & Chiang, 2023). This mechanism can operate even if hot Jupiters form
_in situ_ (Batygin et al., 2016). Such resonant excitations occur more readily
in systems hosting Jupiter-mass planets, and would not operate as efficiently
to create misalignment in systems with smaller planets.

Second, hot Jupiters have multiple possible formation models. While some hot
Jupiters are thought to form either _in situ_ or via dynamically quiet disk
migration (e.g. Becker et al., 2015; Zhu et al., 2018; Hord et al., 2022; Sha
et al., 2023), recent evidence suggests that the majority of hot Jupiters form
via high-eccentricity migration (Zink & Howard, 2023), a process which
requires an initial dynamical excitation due to a companion and subsequent
orbital circularization due to tides raised on the planet’s and star’s
surfaces (Dawson & Johnson, 2018). Many models of hot Jupiter migration invoke
an inclined companion to start the initial migration process, and many of
these mechanisms start with an inclined companion and subsequent Lidov–Kozai
oscillations which perturb the Jupiter into a small-periastron orbit that can
then tidally circularize (Wu & Murray, 2003; Naoz et al., 2012; Petrovich,
2015a; Naoz, 2016; Vick et al., 2019). This picture would produce the lack of
binary companion alignment seen in the systems hosting hot Jupiters in our
sample, as small planets do not have similar origin pathways. However, Lidov-
Kozai oscillations likely explain the origin of only a small fraction (20% or
less; Petrovich, 2015a) of the known hot Jupiters. Instead, Zink & Howard
(2023) suggests that the coplanar high-eccentricity migration mechanism of
Petrovich (2015b) is likely dominant, meaning that even initially coplanar
outer companions can perturb cold Jupiters into hot orbits, where they are
more readily discovered.

Winn et al. (2010) and Rice et al. (2022) find that there is a population of
cool stars hosting hot Jupiters with low stellar obliquities. This population
generally comprises most observed hot Jupiters (Dong & Foreman-Mackey, 2023),
and likely composes a significant fraction of our systems given that 72% of
the hot Jupiters in our sample orbit cool stars. The hot Jupiters with cool
stars probably have predominately aligned stellar spins and planetary orbits,
but we observe a distribution of inclinations consistent with isotropy. These
two observations imply that the distributions of stellar spins and binary
orbits should be roughly uncorrelated.

However, one puzzle remains: under the assumption that we are mostly observing
this population, if we were to measure the obliquity (stellar spin/planet
orbit) angle of the systems in our sample, it should be uncorrelated with the
inclination angle (planet orbit/binary orbit). This is in tension with Behmard
et al. (2022), who found tentative evidence for correlations between obliquity
and inclination distributions. However, Rice et al. (2024) used a larger
sample of systems and did not find evidence for correlations between obliquity
and inclination distributions.

Spin-orbit measurements of the stars hosting hot Jupiters in our sample could
help resolve this uncertainty and determine the angular momentum hierarchies
in these hot-Jupiter hosting systems, enabling a more concrete conclusion on
the dynamical origins of these hot Jupiters.

##  6 Conclusions

Previously, several studies (Christian et al., 2022; Dupuy et al., 2022;
Behmard et al., 2022; Lester et al., 2023) found evidence that the orbits of
exoplanets tend to be aligned with the orbits of wide binary companions in
these systems. We extend the work of Christian et al. (2022) in particular,
significantly increase the size of the sample, and update their conclusions.
We do this by determining the probability distributions of the orbital
parameters of binary systems that host transiting exoplanets (from the TESS
mission) using astrometric measurements from Gaia DR3.

We confirm and amend the findings of these previous papers and show that the
alignment seen previously is dominated by systems containing small planets
(smaller than about 6 $R_{\oplus}$) in binary systems with semimajor axes less
than about 700 AU. The significance of the alignment we detect is considerably
greater than that reported by previous studies.

One interesting result of our work is that at first glance, it appears that
our sample of large planet systems (dominated by hot Jupiters) in wide
binaries are not preferentially aligned to the same degree as small planets.
One must be careful in confidently claiming this to be true given several
selection biases, most importantly the relatively high false positive rate of
hot Jupiters (Section 4.2.2). Despite these hurdles, a bootstrap analysis
shows that the higher false positive rate of hot Jupiters compared to small
planets most likely cannot explain the difference in alignment between the two
samples. This result suggests that the orbits of hot Jupiters in relatively
close ($a\lesssim 700$AU) visual binary systems show less alignment with the
binaries than similar systems containing small planets (Section 4.3.1). This
is likely the signature of the dynamical upheaval believed to be associated
with hot Jupiter migration.

As time goes on, it should be possible to continue strengthening the
conclusions from this work. Future Gaia data releases should result in more
discoveries of binary systems and better constraints on their orbital
parameters. Meanwhile, as TESS continues to observe, and astronomers follow-up
the planet candidates and eliminate more false positives, the purity and size
of the TESS planet sample should improve. We expect that as our sample grows
and improves, the question over the difference in small planet and hot Jupiter
alignment will be resolved with greater certainty.

## Acknowledgements

We acknowledge helpful conversations with Bekki Dawson. S.C. acknowledges
support from the MIT Undergraduate Research Opportunities Program (UROP). The
authors acknowledge the MIT SuperCloud and Lincoln Laboratory Supercomputing
Center for providing high-performance computing resources that have
contributed to the research results reported within this paper. F.J.P
acknowledges financial support from the grant CEX2021-001131-S funded by
MCIN/AEI/ 10.13039/501100011033 and through projects PID2019-109522GB-C52 and
PID2022-137241NB-C43. This paper includes data collected by the TESS mission,
which are publicly available from the Mikulski Archive for Space Telescopes
(MAST). Funding for the TESS mission is provided by NASA’s Science Mission
directorate. This work has made use of data from the European Space Agency
(ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia
Data Processing and Analysis Consortium (DPAC;
https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has
been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement. This research has made use
of NASA’s Astrophysics Data System and the NASA Exoplanet Archive operated by
the California Institute of Technology, under contract with NASA under the
Exoplanet Exploration Program.

## Data Availability

Parameter summaries for binaries and exoplanets in the sample, along with the
entire control sample, are available as supplementary tables. Individual
samples for systems are available from the corresponding author upon request.

## References

  * Albrecht et al. (2012) Albrecht S., et al., 2012, ApJ, 757, 18
  * Albrecht et al. (2022) Albrecht S. H., Dawson R. I., Winn J. N., 2022, arXiv e-prints, p. arXiv:2203.05460
  * Anglada-Escudé et al. (2016) Anglada-Escudé G., et al., 2016, Nature, 536, 437
  * Bailer-Jones et al. (2021) Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Demleitner M., Andrae R., 2021, AJ, 161, 147
  * Batygin (2012) Batygin K., 2012, Nature, 491, 418
  * Batygin et al. (2016) Batygin K., Bodenheimer P. H., Laughlin G. P., 2016, ApJ, 829, 114
  * Becker et al. (2015) Becker J. C., Vanderburg A., Adams F. C., Rappaport S. A., Schwengeler H. M., 2015, ApJ, 812, L18
  * Behmard et al. (2022) Behmard A., Dai F., Howard A. W., 2022, AJ, 163, 160
  * Blunt et al. (2017) Blunt S., et al., 2017, AJ, 153, 229
  * Buchhave et al. (2014) Buchhave L. A., et al., 2014, Nature, 509, 593
  * Christian et al. (2022) Christian S., et al., 2022, AJ, 163, 207
  * Dawson & Johnson (2018) Dawson R. I., Johnson J. A., 2018, ARA&A, 56, 175
  * Dong & Foreman-Mackey (2023) Dong J., Foreman-Mackey D., 2023, AJ, 166, 112
  * Doyle et al. (2011) Doyle L. R., et al., 2011, Science, 333, 1602
  * Dupuy et al. (2022) Dupuy T. J., Kraus A. L., Kratter K. M., Rizzuto A. C., Mann A. W., Huber D., Ireland M. J., 2022, MNRAS, 512, 648
  * El-Badry & Rix (2018) El-Badry K., Rix H.-W., 2018, MNRAS, 480, 4884
  * El-Badry et al. (2021) El-Badry K., Rix H.-W., Heintz T. M., 2021, arXiv e-prints, p. arXiv:2101.05282
  * Epstein-Martin et al. (2022) Epstein-Martin M., Becker J., Batygin K., 2022, ApJ, 931, 42
  * Fabrycky & Tremaine (2007) Fabrycky D., Tremaine S., 2007, ApJ, 669, 1298
  * Fouesneau et al. (2022) Fouesneau M., et al., 2022, arXiv e-prints, p. arXiv:2206.05992
  * Gerbig et al. (2024) Gerbig K., Rice M., et al., 2024, in prep
  * Guerrero et al. (2021) Guerrero N. M., et al., 2021, ApJS, 254, 39
  * Hjorth et al. (2021) Hjorth M., Albrecht S., Hirano T., Winn J. N., Dawson R. I., Zanazzi J. J., Knudstrup E., Sato B., 2021, Proceedings of the National Academy of Science, 118, e2017418118
  * Hord et al. (2022) Hord B. J., et al., 2022, AJ, 164, 13
  * Huang et al. (2016) Huang C., Wu Y., Triaud A. H., 2016, The Astrophysical Journal, 825, 98
  * Jensen et al. (2020) Jensen E. L. N., Akeson R. L., Hersch A. S., 2020, in Five Years After HL Tau: A New Era in Planet Formation. , doi:10.5281/zenodo.4300482
  * Kozai (1962) Kozai Y., 1962, AJ, 67, 591
  * Kraus et al. (2016) Kraus A. L., Ireland M. J., Huber D., Mann A. W., Dupuy T. J., 2016, AJ, 152, 8
  * Kunimoto et al. (2022) Kunimoto M., et al., 2022, ApJS, 259, 33
  * Lai (2014) Lai D., 2014, MNRAS, 440, 3532
  * Lester et al. (2023) Lester K. V., et al., 2023, AJ, 166, 166
  * Lidov (1962) Lidov M. L., 1962, Planet. Space Sci., 9, 719
  * Lindegren (2018) Lindegren L., 2018, Re-normalising the astrometric chi-square in Gaia DR2, GAIA-C3-TN-LU-LL-124, http://www.rssd.esa.int/doc_fetch.php?id=3757412
  * Mann et al. (2019) Mann A. W., et al., 2019, ApJ, 871, 63
  * Moe & Kratter (2021) Moe M., Kratter K. M., 2021, MNRAS, 507, 3593
  * Morton (2015) Morton T. D., 2015, isochrones: Stellar model grid package (ascl:1503.010)
  * Mugrauer (2019) Mugrauer M., 2019, MNRAS, 490, 5088
  * Naoz (2016) Naoz S., 2016, ARA&A, 54, 441
  * Naoz et al. (2012) Naoz S., Farr W. M., Rasio F. A., 2012, ApJ, 754, L36
  * Ngo et al. (2016) Ngo H., et al., 2016, ApJ, 827, 8
  * Offner et al. (2023) Offner S. S. R., Moe M., Kratter K. M., Sadavoy S. I., Jensen E. L. N., Tobin J. J., 2023, in Inutsuka S., Aikawa Y., Muto T., Tomida K., Tamura M., eds, Astronomical Society of the Pacific Conference Series Vol. 534, Protostars and Planets VII. p. 275 (arXiv:2203.10066), doi:10.48550/arXiv.2203.10066
  * Orosz et al. (2012) Orosz J. A., et al., 2012, Science, 337, 1511
  * Pearce et al. (2020) Pearce L. A., Kraus A. L., Dupuy T. J., Mann A. W., Newton E. R., Tofflemire B. M., Vanderburg A., 2020, ApJ, 894, 115
  * Petrovich (2015a) Petrovich C., 2015a, ApJ, 799, 27
  * Petrovich (2015b) Petrovich C., 2015b, ApJ, 805, 75
  * Raghavan et al. (2010) Raghavan D., et al., 2010, ApJS, 190, 1
  * Rice et al. (2022) Rice M., Wang S., Laughlin G., 2022, The Astrophysical Journal Letters, 926, L17
  * Rice et al. (2024) Rice M., Gerbig K., Vanderburg A., 2024, arXiv e-prints, p. arXiv:2401.04173
  * Schlaufman (2010) Schlaufman K. C., 2010, ApJ, 719, 602
  * Schlaufman (2015) Schlaufman K. C., 2015, ApJ, 799, L26
  * Sha et al. (2023) Sha L., et al., 2023, MNRAS, 524, 1113
  * Spalding & Batygin (2014) Spalding C., Batygin K., 2014, ApJ, 790, 42
  * Veras (2016) Veras D., 2016, Royal Society Open Science, 3, 150571
  * Vick et al. (2019) Vick M., Lai D., Anderson K. R., 2019, MNRAS, 484, 5645
  * Winn et al. (2010) Winn J. N., Fabrycky D., Albrecht S., Johnson J. A., 2010, ApJ, 718, L145
  * Wu & Murray (2003) Wu Y., Murray N., 2003, ApJ, 589, 605
  * Zanazzi & Chiang (2023) Zanazzi J. J., Chiang E., 2023, MNRAS,
  * Zanazzi & Lai (2018) Zanazzi J. J., Lai D., 2018, MNRAS, 477, 5207
  * Zhou et al. (2019) Zhou G., et al., 2019, AJ, 158, 141
  * Zhu et al. (2018) Zhu W., Dai F., Masuda K., 2018, Research Notes of the American Astronomical Society, 2, 160
  * Zink & Howard (2023) Zink J., Howard A., 2023, arXiv e-prints, p. arXiv:2310.01567
  * von Zeipel (1910) von Zeipel H., 1910, Astronomische Nachrichten, 183, 345