# Lattice determination of the Batalin-Vilkovisky function
and the strong running interaction

A. C. Aguilar  N. Brito  M. N. Ferreira  J. Papavassiliou  O. Oliveira  P. J.
Silva  University of Campinas - UNICAMP, Institute of Physics Gleb Wataghin,
13083-859 Campinas, São Paulo, Brazil  CFisUC, Department of Physics,
University of Coimbra, Coimbra, 3004-516, Portugal  Centre for Mathematical
Sciences, University of Plymouth, Plymouth, PL4 8AA, United Kingdom
Department of Theoretical Physics and IFIC, University of Valencia and CSIC,
E-46100, Valencia, Spain

###### Abstract

The Batalin-Vilkovisky function is a central component in the modern
formulation of the background field method and the physical applications
derived from it. In the present work we report on novel lattice results for
this particular quantity, obtained by capitalizing on its equality with the
Kugo-Ojima function in the Landau gauge. The results of the lattice simulation
are in very good agreement with the predictions derived from a continuum
analysis based on the corresponding Schwinger-Dyson equations. In addition, we
show that an important relation connecting this function with the ghost
propagator is fulfilled rather accurately. With the aid of these results, we
carry out the first completely lattice-based determination of the process-
independent strong running interaction, employed in a variety of
phenomenological studies.

###### keywords:

Lattice QCD , Background Field Method , Schwinger-Dyson equations , strong
running interaction ,

††journal: Physics Letters B

##  1 Introduction

The background field method (BFM) is a powerful framework that permits the
implementation of the gauge-fixing procedure necessary for quantizing gauge
theories without losing explicit gauge invariance [1, 2, 3, 4, 5, 6, 7, 8, 9,
10]. The formulation of non-Abelian gauge theories within this quantization
scheme affords a plethora of advantages, both from the formal as well as the
practical points of view. Thus, in addition to streamlining a variety of
demonstrations related to renormalization, it allows one to tackle efficiently
longstanding challenges, such as the gauge-invariant truncation of the
Schwinger-Dyson equations (SDEs) [11, 12, 13, 14, 15, 16, 17] or the
definition of physically meaningful renormalization-group invariant (RGI)
quantities, see, e.g., [15].

Intrinsic to the BFM formalism is the duplication of the number of gauge
fields, which are distinguished into “background” and "quantum" types [8, 10];
while the former are not integrated over in the path integral and do not
appear in loops, the latter are identified with the standard gauge fields
known from the conventional formulation, e.g., in the linear covariant gauges
[18].

The correlation (Green’s) functions built out of background fields satisfy
linear Slavnov-Taylor identities (STIs) [19, 20], which are naive
generalizations of tree level relations, without deformations originating from
the ghost-sector of the theory. In addition, and more importantly for our
purposes, they are connected to the corresponding Green’s functions involving
quantum fields by exact relations, known as "background-quantum identities"
(BQIs) [21, 22, 23, 15]. These special identities are derived with the aid of
the Batalin-Vilkovisky (BV) formalism [24, 25], where the original BFM action
is extended through the inclusion of suitable anti-fields and sources.

The common ingredient of this infinite tower of identities is a particular
function, denoted by $G(q)$, which we denominate "BV function" [23]. The
precise field-theoretic definition of this function involves a special
combination of an anti-field and a source; however, a more practical relation
expresses $G(q)$ in terms of the conventional gluon and ghost propagators, and
the ghost-gluon kernel, known from the STI satisfied by the three-gluon vertex
[26, 27, 28, 29, 30, 31, 32, 33].

The BV function appears in a variety of non-perturbative applications, such as
the so-called "block-wise" truncations of the SDE series [34, 35, 36, 37, 38],
and the definition of an effective interaction which is both RGI and process-
independent [39, 40], constituting a common component of every two-to-two on-
shell process. In addition, $G(q)$ is related to the inverse of the ghost
dressing function through an exact relation [41, 42, 43, 40], which enforces
the coincidence of distinct versions of effective charges in the deep infrared
[40]. In that sense, the BV function represents an essential component of the
gauge sector of Yang-Mills theories, participating non-trivially in the
description of the associated dynamics.

Particularly pivotal to the present analysis is the Landau-gauge equality
between the BV function and the so-called Kugo-Ojima (KO) function [44, 41],
whose field-theoretic definition is especially suitable for a lattice
simulation [45, 46, 47, 48]. Exploiting this key relation, we carry out a
large-volume simulation of the KO function, and thus, we simultaneously obtain
the lattice result for the BV function.

The lattice results are contrasted with those obtained from the corresponding
SDE that governs the evolution of $G(q)$, finding rather good agreement. In
addition, the aforementioned relation between $G(q)$ and the ghost dressing
function is shown to be fulfilled at a reasonable level of accuracy.
Furthermore, the lattice results for $G$, together with the data for the gluon
propagator obtained from the same lattice, allow us to present the first
purely lattice-based construction of the effective interaction introduced in
[49, 50].

##  2 The ubiquitous Batalin-Vilkovisky function

Within the BFM approach, the gauge field $A$ appearing in the classical action
is decomposed as $A=B+Q$, where $B$ and $Q$ are the background and quantum
(fluctuating) fields, respectively. Background fields participate in Feynman
diagrams only as external legs, while loops are comprised exclusively by
quantum fields [8, 10]. The key property of the BFM is that the gauge-fixing
may be implemented without compromising explicit gauge invariance.
Specifically, instead of the gauge-fixing term

${\cal L}_{\mathrm{gf}}=\frac{1}{2\xi}(\partial_{\mu}Q^{a\mu})^{2}\,,$ (1)

of the conventional covariant ($R_{\xi}$) gauges [18], one uses

${\widehat{\cal
L}}_{\mathrm{gf}}=\frac{1}{2\xi_{Q}}({\hstretch{2}{\hat{\hstretch{.5}{D}}}}_{\mu}^{ab}Q^{b\mu})^{2}\,,\qquad{\hstretch{2}{\hat{\hstretch{.5}{D}}}}_{\mu}^{ab}=\partial_{\mu}\delta^{ab}+gf^{amb}B_{\mu}^{m}\,,$
(2)

where $f^{abc}$ are the structure constants of SU$(3)$. The key feature that
converts the BFM into a powerful quantization scheme is that the gauge-fixing
choice of Eq. (2) gives rise to a gauge-fixed action which is invariant under
gauge transformations of the background field, $\delta
B^{a}_{\mu}=-g^{-1}\partial_{\mu}\theta^{a}+f^{abc}\theta^{b}B^{c}_{\mu}$,
where $g$ denotes the gauge coupling.

A profound consequence of this invariance is the form of the resulting STIs.
Specifically, when Green’s functions are contracted by the momentum carried by
one of their background gluons, they satisfy ghost-free STIs, akin to the
Takahashi identities [51, 52] known from Abelian theories, such as QED. These
special STIs have far-reaching consequences for renormalization, because they
impose QED-like constraints among the various renormalization constants. In
particular, the wave-function renormalization constant of the $B$ field,
${Z}_{B}$, and the gauge coupling renormalization constant, $Z_{g}$, defined
as

$B^{\,a\mu}_{{\rm{\scriptscriptstyle
R}}}={Z}_{B}^{-1/2}B^{\,a\mu}\,,\qquad\qquad g_{{\rm{\scriptscriptstyle
R}}}=Z_{g}^{-1}g\,,$ (3)

satisfy the key relation (see, e.g., [8, 10])

$Z_{g}={Z}_{B}^{-1/2}\,,$ (4)

in close analogy to the textbook relation $Z_{e}=Z_{A}^{-1/2}$, obeyed by the
renormalization constants of the electric charge and the photon [53].

In what follows we will identify the quantum gauge-fixing parameter $\xi_{Q}$
of the BFM [see Eq. (2)] with the corresponding parameter $\xi$ of the
covariant gauges [see Eq. (1)], i.e., $\xi_{Q}=\xi$; in particular, in the
Landau gauge that we adopt throughout this work, $\xi_{Q}=\xi=0$. With this
identification, the propagator
$\Delta^{ab}_{\mu\nu}(q)=-i\delta^{ab}\Delta_{\mu\nu}(q)$ that connects two
quantum gluons coincides with the standard gluon propagator of the $R_{\xi}$
gauges; in the Landau gauge, $\Delta_{\mu\nu}(q)$ is completely transverse,
i.e.,

$\Delta_{\mu\nu}(q)=\Delta(q){P}_{\mu\nu}(q)\,,\qquad{P}_{\mu\nu}(q):=g_{\mu\nu}-q_{\mu}q_{\nu}/{q^{2}}\,.$
(5)

In addition to $\Delta^{ab}_{\mu\nu}(q)$, in the BFM we have two more gluon
propagators, one connecting the fields $Q^{a}_{\mu}(q)$ with a
$B^{b}_{\nu}(-q)$, and one connecting $B^{a}_{\mu}(q)$ with a
$B^{b}_{\nu}(-q)$; the latter is essential for the ensuing analysis, and will
be denoted by

$\hstretch{2}{\hat{\hstretch{.5}{\Delta}}}^{ab}_{\mu\nu}(q)=-i\delta^{ab}\hstretch{2}{\hat{\hstretch{.5}{\Delta}}}(q)P_{\mu\nu}(q).$
(6)

Figure 1:  Diagrammatic definitions of the ghost-gluon scattering kernel,
$H_{\nu\mu}(r,p,q)$, and the two-point function $\Lambda_{\mu\nu}(q)$, given
in Eq. (8).

It turns out that the gluon propagators $\Delta(q)$ and
$\hstretch{2}{\hat{\hstretch{.5}{\Delta}}}(q)$ are related by the simplest of
all BQIs, namely [21, 23, 42, 54]

$\Delta(q)=[1+G(q)]^{2}\,\hstretch{2}{\hat{\hstretch{.5}{\Delta}}}(q)\,.$ (7)

$G(q)$ is the $g_{\mu\nu}$ component of a certain two-point function,
$\Lambda_{\mu\nu}(q)$, given by (Minkowski space)

$\displaystyle\Lambda_{\mu\nu}(q)$ $\displaystyle:=$ $\displaystyle
ig^{2}C_{\rm A}\int_{k}\Delta^{\rho}_{\mu}(k)D(k+q)H_{\nu\rho}(-q,k+q,-k)$ (8)
$\displaystyle=$ $\displaystyle\underbrace{G(q)}_{\rm
BV\,function}\\!\\!\\!\\!g_{\mu\nu}+\,L(q)\,\frac{q_{\mu}q_{\nu}}{q^{2}}\,,$

where $C_{\mathrm{A}}$ is the Casimir eigenvalue of the adjoint representation
[$N$ for SU$(N)$], $D^{ab}(q)=i\delta^{ab}D(q)$ is the ghost propagator, and
$H_{\nu\mu}(r,p,q)$ denotes the ghost-gluon kernel, see upper panel of Fig. 1.
Note that, formally, $\Lambda_{\mu\nu}(q)$ is a two-point function of a
background source, $\Omega_{\mu}^{a}$, and a gluon anti-field,
$A_{\nu}^{b\star}$, usually denoted in the literature by
${\rm{I}}\\!\Gamma_{\Omega_{\mu}^{a}A_{\nu}^{b\star}}(q)$; the version given
in Eq. (8) [see lower panel of Fig. 1] relies on the fact that
$\Omega_{\mu}^{a}$ and $A_{\nu}^{b\star}$ may be replaced by the Becchi-Rouet-
Stora-Tyutin (BRST) [55, 56] composite operator to which they are coupled.

Next, consider three-point functions (vertices) of the general type
${\mathbf{\Gamma}}_{BXY}$ and ${\mathbf{\Gamma}}_{AXY}$ where $X$ and $Y$
represent general fields (e.g., $X=Y=A$ in the case of a three-gluon vertex,
or $X=\psi$, $Y={\bar{\psi}}$ for the quark-gluon vertex). The vertex BQI has
the general form

${\mathbf{\Gamma}}_{BXY}=[1+G(q)]\,{\mathbf{\Gamma}}_{AXY}+\cdots\,,$ (9)

where the ellipsis denotes terms that involve additional auxiliary functions,
and vanish "on shell".

Higher correlation functions are related by similar, albeit increasingly more
complicated BQIs, with the BV function playing always a prominent role [36].

We next turn to a special identity that relates the two form factors $G(q)$
and $L(q)$ of $\Lambda_{\mu\nu}(q)$ with the ghost dressing function, $F(q)$,
defined as

$F(q)=q^{2}D(q)\,.$ (10)

Specifically, in the Landau gauge, we have the exact relation [40, 54, 49]

$F^{-1}(q)=1+G(q)+L(q)\,.$ (11)

Due to the fact that the function $L(q)$ vanishes at the origin, $L(0)=0$
[40], from Eq. (11) we obtain the important result

$F^{-1}(0)=1+G(0)\,.$ (12)

Note that this last relation enforces the coincidence between the Taylor [57,
12, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67] and "pinch-technique" [68, 39, 15]
effective charges in the deep infrared [69, 70]

##  3 Schwinger-Dyson equations for $G(q)$ and $L(q)$

Within the standard framework of the SDEs, the dynamical equations that govern
the momentum evolution of $G(q)$ and $L(q)$ may be deduced directly from Eq.
(8) by projecting out the $g_{\mu\nu}$ and $q_{\mu}q_{\nu}/q^{2}$ components,
respectively [40]; the tensorial decomposition [32]

$H_{\nu\mu}(r,p,q)=g_{\nu\mu}A_{1}+r_{\mu}r_{\nu}A_{2}+q_{\mu}q_{\nu}A_{3}+r_{\mu}q_{\nu}A_{4}+r_{\nu}q_{\mu}A_{5}\,,$
(13)

with $A_{i}:=A_{i}(r,p,q)$, is also employed.

Setting $f:=1-(q\cdot k)^{2}/(q^{2}k^{2})$ and
${\widetilde{A}}_{i}:=A_{i}(-q,k+q,-k)$, we obtain

$\displaystyle 1+G(q)=$ $\displaystyle\,Z_{c}+\frac{ig^{2}C_{\rm
A}}{3}\int_{k}\Delta(k)D(k+q)K_{G}\,,$ $\displaystyle L(q)=$
$\displaystyle\,\frac{ig^{2}C_{\rm A}}{3}\int_{k}\Delta(k)D(k+q)K_{L}\,,$ (14)

where

$\displaystyle K_{G}=$ $\displaystyle\,(3-f)\,{\widetilde{A}}_{1}-(k\cdot
q)\,f\,{\widetilde{A}}_{4}\,,$ $\displaystyle K_{L}=$
$\displaystyle\,(3-4f)\,{\widetilde{A}}_{1}-f\left[3q^{2}{\widetilde{A}}_{2}+4(q\cdot
k)\,{\widetilde{A}}_{4}\right]\,,$ (15)

and

$\int_{k}:=\frac{1}{(2\pi)^{4}}\\!\int_{-\infty}^{+\infty}\\!\\!\mathrm{d}^{4}k$
(16)

is the properly regularized measure.

Note that, due to the validity of Eq. (11), the $G(q)$ in Eq. (14) is
renormalized through the $Z_{c}$, namely the ghost renormalization constant.
For its determination, one must consider the corresponding SDE for the ghost
dressing function $F(q)$, see, e.g., [38]. This dynamical equation involves
the ghost-gluon vertex as one of its ingredients; this latter vertex is finite
in the Landau gauge, and its attendant renormalization constant, $Z_{1}$, is
independent of the cutoff. Within the Taylor scheme [62, 71, 72] that we will
employ, we have $Z_{1}=1$, and the SDE for $F(q)$ assumes the form

$F^{-1}(q)=Z_{c}+ig^{2}C_{\rm
A}\int_{k}f\Delta(k)D(k+q){\widetilde{B}}_{1}\,,$ (17)

where ${\widetilde{B}}_{1}:=B_{1}(-q,k+q,-k)$, with $B_{1}$ denoting the
classical form factor of the ghost-gluon vertex. Then, the expression for
$Z_{c}$ to be used in Eq. (14) is obtained from Eq. (17) by imposing the
renormalization condition $F(\mu)=1$, where $\mu$ is the renormalization
point; throughout this work we use $\mu=4.3$ GeV.

The numerical treatment of the SDE system formed by Eqs. (14) and (17)
proceeds along the lines described in [73, 40], using up-to-date ingredients
as external inputs. Specifically, we employ general kinematics results for
${\widetilde{B}}_{1}$ and ${\widetilde{A}}_{i}$, determined from their own
SDEs in [32]; instead, in [73, 40] these form factors were set at their tree-
level values. Moreover, we employ directly a fit to the lattice gluon
propagator $\Delta(q)$ (lower panel of Fig. 4), obtained from the same lattice
setups used for the determination of $1+G(q)$ and $F(q)$. The resulting SDE
solutions for $G(q)$ and $L(q)$ will be discussed in Sec. 5.

##  4 Lattice simulation: theory and setup

Both the formal definition of $G(q)$ in terms of BV fields and the alternative
presented in Eq. (8) are unsuitable for performing lattice simulations.
Instead, we will take advantage of the known equality between $G(q)$ and the
KO function111 This function is associated with a standard confinement
criterion [44, 41]; for a review, see [12] . It is well-known that the $u$
simulated on the lattice does not comply with this criterion [45, 46], nor do
the ”decoupling” propagators, see, e.g., [74, 75, 76, 77] and [78, 37, 79].,
to be denoted by $u(q)$; specifically, in the Landau gauge, we have [41, 42,
43, 73].

$G(q)=u(q)\,.$ (18)

The field-theoretic definition of $u(q)$ that has been implemented on the
lattice involves two composite operators,

${\cal A}^{a}_{\mu}(x):=D_{\mu}^{ae}(x)\,c^{e}(x)\,,\qquad{\cal
B}^{a}_{\nu}(x):=f^{bcd}A_{\nu}^{d}(x)\,\bar{c}^{c}(x)\,,$ (19)

where $D_{\mu}^{ab}(x)=\partial_{\mu}\delta^{ac}+gf^{amb}A^{m}_{\mu}(x)$ is
the covariant derivative in the adjoint representation, and the two-point
function $\mathcal{U}^{ab}_{\mu\nu}(q)$, defined as (Euclidean space)

$\mathcal{U}^{ab}_{\mu\nu}(q)=\int d^{4}xe^{iq(x-y)}\ \langle 0|T\left({\cal
A}^{a}_{\mu}(x)\,{\cal B}^{b}_{\nu}(y)\right)|0\rangle\,,$ (20)

where $T$ denotes the standard time-ordering operation. The KO function is the
scalar co-factor of $\mathcal{U}^{ab}_{\mu\nu}(q)$,

$\mathcal{U}^{ab}_{\mu\nu}(q)=\delta^{ab}\left(\delta^{\mu\nu}-\frac{q_{\mu}q_{\nu}}{q}\right)u(q)\,.$
(21)

An explicit lattice definition for $\mathcal{U}^{ab}_{\mu\nu}(q)$ is given by

$\mathcal{U}^{ab}_{\mu\nu}(q)\\!=\\!\frac{1}{V}\left\langle\sum_{x,y,z}\sum_{c,d,e}e^{-iq\cdot(x-y)}\left(D_{\mu}\right)^{ae}\\!\\!\\!(x;z)\left(M^{-1}\right)^{ec}\\!\\!\\!(z;y)f^{bcd}A_{\nu}^{d}(y)\right\rangle_{U}$
(22)

where $\left\langle\ldots\right\rangle_{U}$ denotes the Monte-Carlo sample
average over the gauge field configurations $U$, and $M$ is a suitably
discretized version of the Faddeev-Popov operator,
$\partial_{\mu}D^{ab}_{\mu}$; note that the inverse of $M$ is related to the
ghost propagator by

$D^{ab}(q)=\int d^{4}xe^{iq(x-y)}\ \langle
0|T\left(M^{-1}\right)^{ab}\\!\\!\\!(x;y)|0\rangle\,.$ (23)

With the aid of Eq. (21), the scalar function $u(q)$ is given by

$u(q)=\frac{1}{24}\sum_{\mu,a}\mathcal{U}^{aa}_{\mu\mu}(q)\,.$ (24)

In order to study the KO function on the lattice, we rely on Eq. (22).
However, for practical reasons, it is convenient to compute
$\mathcal{U}^{ab}_{\mu\nu}(q)$ using a point source $y_{0}$ in the inversion
of the lattice Faddeev-Popov operator, $M$,

$\mathcal{U}^{ab}_{\mu\nu}(q\\!=\\!\left\langle\sum_{x,z}\sum_{c,d,e}e^{-iq\cdot(x-y_{0})}\left(D_{\mu}\right)^{ae}\\!\\!\\!(x;z)\left(M^{-1}\right)^{ec}\\!\\!\\!(z;y_{0})f^{bcd}A_{\nu}^{d}(y_{0})\right\rangle_{U}$
(25)

where the lattice definitions for $D$, $M$, and $A$ can be found in [46, 80,
81].

The most important points of the lattice simulation may be summarized as
follows:

($\it i$) The gauge field, $A_{\mu}^{a}$, covariant derivative,
$D_{\mu}^{ab}$, and Faddeev-Popov operator, $M_{xy}^{ab}$, are suitably
discretized according to [46, 47, 77, 82].

($\it ii$) The Landau gauge is fixed by identifying the configurations
belonging to the Gribov region; this is achieved through the maximization of a
certain functional, whose discretized form is given in [46, 47]. Numerically,
the required maximization is performed with the Fourier Accelerated Steepest
Descent method [83].

($\it iii$) Since $M_{xy}^{ab}$ is a singular operator, it can only be
inverted in the subspace orthogonal to its zero modes. Nevertheless, this is
sufficient for computing $u(q)$ for $q\neq 0$ [81]. In this case, the
inversion problem in the orthogonal subspace can be transformed into a sparse
linear system of equations, as explained in detail in [46, 47, 81]. Then,
since $M_{xy}^{ab}$ is real and symmetric, the resulting system can be solved
efficiently with the conjugate gradient method [84].

($\it iv$) We consider quenched lattice ensembles generated with the Wilson
gauge action [85], with $\beta=6.0$. In this case, the Sommer parameter method
[86] yields a lattice spacing of $a=0.0962$ fm ($a^{-1}=2.05$ GeV) [87]. We
employ two setups, with lattice volumes $V=64^{4}$, and $V=80^{4}$, in lattice
units; the corresponding physical volumes extend from $~{}(6\textrm{ fm})^{4}$
to $~{}(8\textrm{ fm})^{4}$. Further details on sampling, gauge fixing, and
definitions can be found in [88].

($\it v$) The number of gauge field configurations considered for the present
study was $700$ for $V=64^{4}$, and $500$ for $V=80^{4}$.

($\it vi$) The computer simulations were performed with the help of the Chroma
library [89], using PFFT [90] for the necessary Fourier transforms.

##  5 Results

Even though the lattice simulation uses the expression in Eq. (21) as its
starting point, in what follows we present the results in terms of $G(q)$,
which is the focal point of this study, rather than $u(q)$. Evidently, by
virtue of Eq. (18), the results may be recast trivially in terms of $u(q)$, by
setting $G\\!\to\\!u$ throughout.

Let us denote the bare $G(q)$ obtained on the lattice by
$G_{{\rm{\scriptscriptstyle B}}}(q,a)$, making manifest the dependence on the
lattice spacing $a$. The result for $G_{{\rm{\scriptscriptstyle B}}}(q,a)$
from our simulations are shown in Fig. 2, for both lattice volumes. We note
that the volume dependence is rather weak; indeed, the data from both volumes
agree within errors.

Ideally, in order to obtain a reliable comparison with the SDE results, the
$G_{{\rm{\scriptscriptstyle B}}}(q,a)$ should be extrapolated to the continuum
limit, $a\to 0$, and be renormalized using the Taylor scheme. Instead, in what
follows we adopt a slightly modified version of the procedure developed in
[46], whereby $G_{{\rm{\scriptscriptstyle B}}}(q,a)$ is renormalized by
imposing Eq. (11).

Figure 2: Bare lattice results for the BV function,
$G_{{\rm{\scriptscriptstyle B}}}(q,a)$.

To begin with, the continuum limit $G_{{\rm{\scriptscriptstyle
B}}}(q):=G_{{\rm{\scriptscriptstyle B}}}(q,0)$ is related to the lattice
$G_{{\rm{\scriptscriptstyle B}}}(q,a)$ by

$G_{{\rm{\scriptscriptstyle B}}}(q,a)=G_{{\rm{\scriptscriptstyle
B}}}(q)+s(q)\,,$ (26)

where $s(q)$ accounts for lattice artifacts, and, in principle, depends on the
momentum $q$. However, the momentum dependence of $s(q)$ cannot be determined
by the analysis employed below. Therefore, in order to proceed, we assume that
the dependence of $s(q)$ on the momentum is fairly mild [77] and treat it as a
constant, i.e., $s(q)\to s$.

Then, the renormalized lattice $1+G(q)$ in the Taylor scheme is given by

$1+G(q)=Z_{c}[1+G_{{\rm{\scriptscriptstyle B}}}(q,a)-s]\,,$ (27)

where the parameter $s$ acts as an effective subtractive renormalization
constant. We emphasize that subtractive terms in the renormalization have been
shown to improve the agreement with the continuum theory in studies of the
quark propagator and the quark-gluon vertex [91, 92].

The value of $Z_{c}$ can be determined formally by imposing a renormalization
condition, i.e., a value for the renormalized $1+G(\mu)$, at some scale $\mu$.
However, since in the Taylor scheme $F(\mu)=1$, and $L(\mu)$ generally does
not vanish, Eq. (11) implies that we cannot impose $1+G(\mu)=1$.

Instead, we determine the value of $1+G(\mu)$ by first computing $1+G(q)$ at
one loop,

$1+G^{(1)}(q)=1+\frac{\alpha_{s}C_{\rm
A}}{16\pi}\left[3\ln(q^{2}/\mu^{2})-2\right]\,,$ (28)

with $\alpha_{s}=g^{2}/4\pi$. Then, using the value of $\alpha_{s}=0.216$ for
the Taylor coupling at $\mu=4.3$ GeV [62] we obtain

$1+G^{(1)}(\mu)=0.974\,.$ (29)

Interestingly, the SDE result yields $1+G(\mu)=0.973$, indicating that the
perturbative regime of $1+G(q)$ has been safely reached for $q=\mu=4.3$ GeV.
Hence, we can fix the scale of the lattice $1+G(q)$ by imposing
$1+G(\mu)=1+G^{(1)}(\mu)$, such that

$Z_{c}=\frac{1+G^{(1)}(\mu)}{1+G_{{\rm{\scriptscriptstyle
B}}}(\mu,a)-s}=\frac{0.974}{0.922-s}\,,$ (30)

where used $G_{{\rm{\scriptscriptstyle B}}}(\mu,a)=-0.078$, obtained from the
bare data of Fig. 2.

At this point, $Z_{c}$ still depends on the unknown constant $s$. To fix its
value, we impose the validity of Eq. (11), i.e.,

$Z_{c}\left[1+G_{{\rm{\scriptscriptstyle
B}}}(q,a)-s\right]+L(q)=F_{\\!\\!{\rm{\scriptscriptstyle lat}}}^{-1}(q)\,,$
(31)

where $F_{\\!\\!{\rm{\scriptscriptstyle lat}}}(q)$ is the renormalized lattice
ghost dressing function, obtained self-consistently from the same lattice
setups used to compute $G_{{\rm{\scriptscriptstyle B}}}(q,a)$.

Now, the function $L(q)$ is not directly computed on the lattice.
Nevertheless, for small $q$ we can assume for $L(q)$ a Taylor expansion,

$L(q)=a_{1}q^{2}+a_{2}q^{4}+{\cal O}(q^{6})\,,$ (32)

whose $0$-th order term vanishes on account of $L(0)=0$.

Then, we can determine the parameters $s$, $a_{1}$, and $a_{2}$, by rephrasing
Eq. (31) as a $\chi^{2}$ minimization problem. Specifically, we minimize the
$\chi^{2}$ defined by

$\chi^{2}:=\sum_{i}\left\\{Z_{c}\left[1+G_{{\rm{\scriptscriptstyle
B}}}(q_{i},a)-s\right]+a_{1}q_{i}^{2}+a_{2}q_{i}^{4}-F_{\\!\\!{\rm{\scriptscriptstyle
lat}}}^{-1}(q_{i})\right\\}^{2}\,,$ (33)

where $q_{i}$ are the lattice points, and $Z_{c}$ is given by Eq. (30).

Figure 3:  Top: Results for $1+G(q)$ obtained from the lattice (points) and
the SDE prediction (red curve). Bottom: Comparison of $[1+G(q)]^{-1}$
(circles/squares) and $F(q)$ (stars/triangles).

Since Eq. (33) is the result of a Taylor expansion, it is supposed to be used
only for small $q$; specifically, we choose the points $q_{i}<0.5$ GeV.
Consequently, there are too few points to reliably determine the parameters
$s$, $a_{1}$ and $a_{2}$ for the two lattice setups individually. Instead, we
apply the above procedure to the _combined_ data points of both volumes, for
which there are $18$ points in the fitting window. We have verified that for
the combined data, the value of $s$ is stable against varying the order of the
Taylor expansion in Eq. (32), as well as the fitting window. Specifically, we
obtain $s=0.121$, $a_{1}=0.458~{}\text{GeV}^{-2}$, and
$a_{2}=-0.905~{}\text{GeV}^{-4}$, for which $\chi^{2}=2\times 10^{-3}$.

Using the above value of $s$ in Eqs. (27) and (30), we obtain the renormalized
lattice $1+G(q)$, shown as points on the top panel of Fig. 3. On the same
panel we show also the SDE prediction (red continuous), finding an excellent
agreement. On the bottom panel of Fig. 3, we compare $F(q)$ to
$[1+G(q)]^{-1}$; evidently, Eq. (12) is satisfied within errors, as a
consequence of the renormalization procedure employed.

The comparison between lattice and SDE results may be taken a step further, by
considering the full momentum dependence of the form factor $L(q)$. In
particular, a lattice-derived result for $L(q)$ may be obtained from Eq. (11),
by substituting in it the lattice data for $F(q)$ and $1+G(q)$. The outcome of
this procedure is shown in Fig. 4; the comparison with the SDE prediction (red
continuous curve) reveals a reasonable agreement. The blue curve represents a
fit to the lattice $L(q)$, given by the simple functional form

$L(q)=\frac{q^{2}/c_{1}^{2}+(q^{2}/c_{2}^{2})^{2}}{1+q^{2}/t_{1}^{2}+(q^{2}/t_{2}^{2})^{2}\ln^{d_{{\rm{\scriptscriptstyle
L}}}}(1+q^{2}/\Lambda_{{\rm{\scriptscriptstyle T}}}^{2})}\,,$ (34)

where $\Lambda_{{\rm{\scriptscriptstyle T}}}=425$ MeV is the value of
$\Lambda_{{\rm{\scriptscriptstyle QCD}}}$ in the Taylor scheme [62], and
$d_{{\rm{\scriptscriptstyle L}}}=35/44$ is the anomalous dimension of $L(q)$
[93]. The fitting parameters are given by $c_{1}^{2}=2.18~{}\text{GeV}^{2}$,
$c_{2}^{2}=0.575~{}\text{GeV}^{2}$, $t_{1}^{2}=0.0981~{}\text{GeV}^{2}$, and
$t_{2}^{2}=0.188~{}\text{GeV}^{2}$. Note that Eq. (34) satisfies $L(0)=0$ by
construction.

Figure 4:  Upper panel: $L(q)$ obtained from the lattice $1+G(q)$ and $F(q)$
of Fig. 3 (points) compared to the SDE prediction (red continuous). Also shown
is a fit given by Eq. (34) (blue dashed). Lower panel: The gluon propagator
simulated on the same lattices as the BV function [94].

##  6 RGI running interaction strength

As has been explained in a series of articles [15, 40, 49], the fundamental
relation of Eq. (4) allows for the definition of a propagator-like RGI
quantity, exactly as happens in QED when the photon vacuum polarization is
multiplied by $e^{2}$. In particular, since
$\hstretch{2}{\hat{\hstretch{.5}{\Delta}}}(q)$ renormalizes through the
$Z_{B}$ introduced in Eq. (4), i.e.,

$\hstretch{2}{\hat{\hstretch{.5}{\Delta}}}_{{\rm{\scriptscriptstyle
R}}}(q)=Z_{B}^{-1}\hstretch{2}{\hat{\hstretch{.5}{\Delta}}}(q)\,,$ (35)

the combination

${\widehat{d}}(q):=\alpha_{s}{\hstretch{2}{\hat{\hstretch{.5}{\Delta}}}(q)}\,,$
(36)

is RGI by virtue of Eq. (4). Indeed, ${\widehat{d}}(q)$ retains exactly the
same form before and after renormalization, and, consequently, does not depend
on the renormalization point $\mu$, nor on the renormalization scheme
employed.

The BV function enters into the definition of ${\widehat{d}}(q)$ when the
central relation in Eq. (7) is invoked,

${\widehat{d}}(q):=\frac{\alpha_{s}\Delta(q)}{[1+G(q)]^{2}}\,;$ (37)

in this form, ${\widehat{d}}(q)$ is known in the literature as the “RGI
running interaction strength” [49].

Figure 5: Diagrammatic representation of the BQI-induced rearrangement in the
case of a generic $X_{1}X_{2}\to X_{3}X_{4}$ scattering, corresponding to the
S-matrix element ${\cal T}_{X_{1}X_{2}\to X_{3}X_{4}}$ of Eq. (38).

What is particularly interesting about ${\widehat{d}}(q)$ is that it does not
represent a mere field-theoretic construct, but admits a clear physical
interpretation. Specifically, ${\widehat{d}}(q)$ constitutes a common
component of any two-to-two physical processes, $X_{1}(k_{1})X_{2}(k_{2})\to
X_{3}(k_{3})X_{4}(k_{4})$, where an off-shell gluon, carrying momentum
$q=k_{1}-k_{2}=k_{3}-k_{4}$, is exchanged, see Fig. 5. In particular, with the
aid of the BQIs in Eqs. (7) and (9), the S-matrix element ${\cal
T}_{X_{1}X_{2}\to X_{3}X_{4}}$ may be cast in the form

$\displaystyle{\cal T}_{X_{1}X_{2}\to X_{3}X_{4}}=$
$\displaystyle\,\,\left\\{g\,{\mathbf{\Gamma}}_{\\!AX_{1}X_{2}}\right\\}\Delta(q)\left\\{g\,{\mathbf{\Gamma}}_{AX_{3}X_{4}}\right\\}$
$\displaystyle=$
$\displaystyle\,\,\left\\{g\,[1+G(q)]^{-1}{{\mathbf{\Gamma}}}_{\\!BX_{1}X_{2}}\right\\}\Delta(q)\left\\{g\,[1+G(q)]^{-1}{{\mathbf{\Gamma}}}_{\\!BX_{3}X_{4}}\right\\}$
$\displaystyle=$
$\displaystyle\,\,{{\mathbf{\Gamma}}}_{\\!BX_{1}X_{2}}\underbrace{\left\\{g^{2}[1+G(q)]^{-2}\Delta(q)\right\\}}_{4\pi\,{\widehat{d}}(q)}{{\mathbf{\Gamma}}}_{\\!BX_{3}X_{4}}\,.$
(38)

Due to the general validity of the BQIs, the steps leading to the last line of
Eq. (38) may be followed regardless of the detailed nature of the initial and
final states. In that sense, ${\widehat{d}}(q)$ captures the process-
independent contribution, common to all such processes, whilst the process-
dependence, i.e., the part that carries the information about the specific
nature of the initial and final states, is encoded in the vertices
${{\mathbf{\Gamma}}}_{\\!BX_{1}X_{2}}$ and
${{\mathbf{\Gamma}}}_{\\!BX_{3}X_{4}}$.

The quantity ${\widehat{d}}(q)$ has mass dimension of $-2$; from it one may
define the dimensionless RGI interaction ${\cal I}(q)$ [49]

${\cal I}(q):=q^{2}{\widehat{d}}(q)\,.$ (39)

As explained in [49], this quantity provides the strength required in order to
describe ground-state hadron observables using SDEs in the matter sector of
the theory. As was argued there, the physics encoded in ${\cal I}(q^{2})$
reconciles nonperturbative continuum QCD with ab initio predictions of basic
hadron properties.

Figure 6:  Lattice results for ${\widehat{d}}(q)$ (top) and ${\cal I}(q)$
(bottom), represented by points. On both panels, the present SDE result (red
continuous curve) and that of [49] (green dotted) are shown for comparison.
The fit to the lattice data, given by Eq. (40) is shown in blue dashed.

In Fig. 6, the red continuous curves correspond to the ${\widehat{d}}(q)$
(top) and ${\cal I}(q)$ (bottom), obtained by combining the $1+G(q)$ from Fig.
3 with the gluon propagator, self-consistently simulated on the same lattice
setups (see lower panel of Fig. 4). For the Taylor coupling at $\mu=4.3$ GeV,
we use $\alpha_{s}=0.216$ [62]; note that, due to the RGI nature of
${\widehat{d}}(q)$, any other set $\\{\mu,\alpha_{s}(\mu)\\}$ yields precisely
the same answer. In addition, the ${\widehat{d}}(q)$ and ${\cal I}(q)$
obtained from the SDE analysis of [49] are displayed as green dotted curves.
Note that the present treatment leads to closer agreement between lattice and
SDE results; we have confirmed that this is mainly due to the non-perturbative
dressing of the form factors ${\widetilde{B}}_{1}$ and ${\widetilde{A}}_{i}$,
discussed in the last paragraph of Sec. 3.

Moreover, in Fig. 6 we show as a blue dashed curve a fit to the lattice
${\widehat{d}}(q)$, given by the functional form

$\displaystyle{\widehat{d}}(q)=$
$\displaystyle\,\frac{{\widehat{d}}(0)\left[1+a_{1}q^{2}+a_{2}q^{4}\right]}{1+b_{1}q^{2}+b_{2}q^{4}+\beta_{0}{\widehat{d}}(0)a_{2}q^{6}\ln(1+q^{2}/\Lambda_{{\rm{\scriptscriptstyle
T}}}^{2})}\,,$ (40)

where $\beta_{0}=11/(4\pi)$ is the first coefficient of the QCD beta function,
$\Lambda_{{\rm{\scriptscriptstyle T}}}=425$ MeV [62], and the fitting
parameters are ${\widehat{d}}(0)=17.1~{}\text{GeV}^{-2}$,
$a_{1}=0.0585~{}\text{GeV}^{-2}$, $a_{2}=0.0274~{}\text{GeV}^{-4}$,
$b_{1}=3.51~{}\text{GeV}^{-2}$, $b_{2}=9.33~{}\text{GeV}^{-4}$. Note that Eq.
(40) enforces the one loop-running of the QCD coupling in the ultraviolet,
i.e., ${\cal I}(q)\to[\beta_{0}\ln(q^{2}/\Lambda_{{\rm{\scriptscriptstyle
T}}}^{2})]^{-1}$ at large $q$.

##  7 Conclusions

The quantitative agreement between the BV function simulated on the lattice
and the corresponding SDE results exposes once again the underlying
consistency of a large number of concepts and techniques, developed over a
period of several years, see, e.g., [95]. Especially interesting in that
regard is the interplay between the ghost and gauge sectors of the theory,
which are nontrivially intertwined by the SDEs. Note, in particular, that the
lattice gluon propagator in Fig. 4 is used as input in the SDEs of Eqs. (14)
and (17), which produce the prediction for $G(q)$ that is subsequently
compared with the corresponding lattice result in Fig. 3. Let us finally
emphasize that the effective interactions ${\widehat{d}}(q)$ and ${\cal I}(q)$
shown in Fig. 6 correspond to the pure Yang-Mills case; for phenomenological
applications they must be modified to include effects from dynamical quarks,
in the spirit of [93, 50, 96]. In fact, the results of the present work
bolster up our confidence in the SDE derivations that lead to the
"unquenching" of these RGI quantities.

## Acknowledgements

The authors acknowledge the computing time provided by the Laboratory for
Advanced Computing at the University of Coimbra (FCT contracts 2021.09759.CPCA
and 2022.15892.CPCA.A2). Work supported by FCT contracts CERN/FIS-
PAR/0023/2021, UIDB/04564/2020, and UIDP/04564/2020. P. J. S. acknowledges
financial support from FCT under Contract No. CEECIND/00488/2017. A. C. A is
supported by the CNPq grant 310763/2023-1, and acknowledges financial support
from project 464898/2014-5 (INCT-FNA). M. N. F. and J. P. are supported by the
Spanish MICINN grant PID2020-113334GB-I00. M. N. F. acknowledges financial
support from Generalitat Valenciana through contract CIAPOS/2021/74. J. P.
also acknowledges funding from the Generalitat Valenciana grant
CIPROM/2022/66. This work was granted access to the HPC resources of the PDC
Center for High Performance Computing at the KTH Royal Institute of
Technology, Sweden, made available within the Distributed European Computing
Initiative by the PRACE-2IP, receiving funding from the European Community’s
Seventh Framework Programme (FP7/ 2007-2013) under grant agreement No.
RI-283493. The use of Lindgren has been provided under DECI-9 project
COIMBRALATT. We acknowledge that the results of this research have been
achieved using the PRACE-3IP project (FP7 RI- 312763) resource Sisu based in
Finland at CSC. The use of Sisu has been provided under DECI-12 project
COIMBRALATT2.

## References

  * DeWitt [1967] B. S. DeWitt, Phys. Rev. 162 (1967) 1195–1239.
  * Honerkamp [1972] J. Honerkamp, Nucl. Phys. B 48 (1972) 269–287.
  * Kallosh [1974] R. E. Kallosh, Nucl. Phys. B 78 (1974) 293–312.
  * Kluberg-Stern and Zuber [1975] H. Kluberg-Stern, J. B. Zuber, Phys. Rev. D 12 (1975) 482–488.
  * Arefeva et al. [1975] I. Y. Arefeva, L. D. Faddeev, A. A. Slavnov, Theor. Math. Phys. 21 (1975) 1165.
  * Abbott [1981] L. Abbott, Nucl. Phys. B 185 (1981) 189–203.
  * Weinberg [1980] S. Weinberg, Phys. Lett. B91 (1980) 51.
  * Abbott [1982] L. F. Abbott, Acta Phys. Polon. B13 (1982) 33.
  * Shore [1981] G. M. Shore, Annals Phys. 137 (1981) 262.
  * Abbott et al. [1983] L. F. Abbott, M. T. Grisaru, R. K. Schaefer, Nucl. Phys. B 229 (1983) 372–380.
  * Roberts and Williams [1994] C. D. Roberts, A. G. Williams, Prog. Part. Nucl. Phys. 33 (1994) 477–575.
  * Alkofer and von Smekal [2001] R. Alkofer, L. von Smekal, Phys. Rept. 353 (2001) 281\.
  * Fischer [2006] C. S. Fischer, J. Phys. G 32 (2006) R253–R291.
  * Roberts [2008] C. D. Roberts, Prog. Part. Nucl. Phys. 61 (2008) 50–65.
  * Binosi and Papavassiliou [2009] D. Binosi, J. Papavassiliou, Phys. Rept. 479 (2009) 1–152.
  * Cloet and Roberts [2014] I. C. Cloet, C. D. Roberts, Prog. Part. Nucl. Phys. 77 (2014) 1–69.
  * Huber [2020] M. Q. Huber, Phys. Rept. 879 (2020) 1–92.
  * Fujikawa et al. [1972] K. Fujikawa, B. W. Lee, A. I. Sanda, Phys. Rev. D6 (1972) 2923–2943.
  * Taylor [1971] J. Taylor, Nucl. Phys. B 33 (1971) 436–444.
  * Slavnov [1972] A. Slavnov, Theor. Math. Phys. 10 (1972) 99–107.
  * Grassi et al. [2001a] P. A. Grassi, T. Hurth, M. Steinhauser, Annals Phys. 288 (2001a) 197–248.
  * Grassi et al. [2001b] P. A. Grassi, T. Hurth, M. Steinhauser, Nucl. Phys. B610 (2001b) 215–250.
  * Binosi and Papavassiliou [2002] D. Binosi, J. Papavassiliou, Phys. Rev. D66 (2002) 025024\.
  * Batalin and Vilkovisky [1977] I. A. Batalin, G. A. Vilkovisky, Phys. Lett. B69 (1977) 309–312.
  * Batalin and Vilkovisky [1983] I. A. Batalin, G. A. Vilkovisky, Phys. Rev. D 28 (1983) 2567–2582. [Erratum: Phys.Rev.D 30, 508 (1984)].
  * Marciano and Pagels [1978] W. J. Marciano, H. Pagels, Phys. Rept. 36 (1978) 137\.
  * Ball and Chiu [1980] J. S. Ball, T.-W. Chiu, Phys. Rev. D 22 (1980) 2550\. [Erratum: Phys. Rev. D 23, 3085 (1981)].
  * Pascual and Tarrach [1984] P. Pascual, R. Tarrach, Lect. Notes Phys. 194 (1984) 1–277.
  * Davydychev et al. [1996] A. I. Davydychev, P. Osland, O. Tarasov, Phys. Rev. D 54 (1996) 4087–4113. [Erratum: Phys. Rev. D 59, 109901 (1999)].
  * Gracey [2011] J. A. Gracey, Phys. Rev. D84 (2011) 085011.
  * Gracey [2014] J. Gracey, Phys. Rev. D 90 (2014) 025014.
  * Aguilar et al. [2019a] A. C. Aguilar, M. N. Ferreira, C. T. Figueiredo, J. Papavassiliou, Phys. Rev. D99 (2019a) 034026.
  * Aguilar et al. [2019b] A. C. Aguilar, M. N. Ferreira, C. T. Figueiredo, J. Papavassiliou, Phys. Rev. D99 (2019b) 094010.
  * Aguilar and Papavassiliou [2006] A. C. Aguilar, J. Papavassiliou, J. High Energy Phys. 12 (2006) 012.
  * Binosi and Papavassiliou [2008a] D. Binosi, J. Papavassiliou, Phys. Rev. D77 (2008a) 061702.
  * Binosi and Papavassiliou [2008b] D. Binosi, J. Papavassiliou, J. High Energy Phys. 11 (2008b) 063.
  * Aguilar et al. [2008] A. C. Aguilar, D. Binosi, J. Papavassiliou, Phys. Rev. D78 (2008) 025010.
  * Aguilar et al. [2023] A. C. Aguilar, M. N. Ferreira, D. Ibañez, B. M. Oliveira, J. Papavassiliou, Eur. Phys. J. C 83 (2023) 86.
  * Binosi and Papavassiliou [2003] D. Binosi, J. Papavassiliou, Nucl. Phys. Proc. Suppl. 121 (2003) 281–284.
  * Aguilar et al. [2009] A. C. Aguilar, D. Binosi, J. Papavassiliou, J. Rodriguez-Quintero, Phys. Rev. D80 (2009) 085018\.
  * Kugo [1995] T. Kugo, in: International Symposium on BRS Symmetry on the Occasion of Its 20th Anniversary, arXiv:hep-th/9511033.
  * Grassi et al. [2004] P. A. Grassi, T. Hurth, A. Quadri, Phys. Rev. D70 (2004) 105014.
  * Kondo [2009] K.-I. Kondo, Phys. Lett. B 678 (2009) 322–330.
  * Kugo and Ojima [1979] T. Kugo, I. Ojima, Prog. Theor. Phys. Suppl. 66 (1979) 1–130.
  * Nakajima and Furui [2000] H. Nakajima, S. Furui, Nucl. Phys. B Proc. Suppl. 83 (2000) 521–523.
  * Sternbeck [2006] A. Sternbeck, The Infrared behavior of lattice QCD Green’s functions, Ph.D. thesis, Humboldt-University Berlin, arXiv:0609016.
  * Brito [2023] N. M. R. Brito, Lattice computation of the Kugo-Ojima correlation function, Master thesis, University of Coimbra, arXiv:2310.09010.
  * Brito et al. [2024] N. Brito, O. Oliveira, P. J. Silva, J. Papavassiliou, M. N. Ferreira, A. C. Aguilar, PoS LATTICE2023 (2024) 380\.
  * Binosi et al. [2015] D. Binosi, L. Chang, J. Papavassiliou, C. D. Roberts, Phys. Lett. B742 (2015) 183–188.
  * Cui et al. [2020] Z.-F. Cui, J.-L. Zhang, D. Binosi, F. de Soto, C. Mezrag, J. Papavassiliou, C. D. Roberts, J. Rodríguez-Quintero, J. Segovia, S. Zafeiropoulos, Chin. Phys. C 44 (2020) 083102\.
  * Ward [1950] J. C. Ward, Phys. Rev. 78 (1950) 182.
  * Takahashi [1957] Y. Takahashi, Nuovo Cim. 6 (1957) 371.
  * Itzykson and Zuber [1980] C. Itzykson, J. B. Zuber, Quantum Field Theory, International Series in Pure and Applied Physics, New York, USA: Mcgraw-Hill (1980) 705 p., 1980.
  * Binosi and Quadri [2013] D. Binosi, A. Quadri, Phys. Rev. D88 (2013) 085036\.
  * Becchi et al. [1976] C. Becchi, A. Rouet, R. Stora, Annals Phys. 98 (1976) 287–321.
  * Tyutin [1975] I. V. Tyutin, LEBEDEV-75-39 (1975).
  * von Smekal et al. [1997] L. von Smekal, R. Alkofer, A. Hauck, Phys. Rev. Lett. 79 (1997) 3591–3594.
  * Fischer and Alkofer [2003] C. S. Fischer, R. Alkofer, Phys. Rev. D67 (2003) 094020\.
  * Bloch et al. [2004] J. C. R. Bloch, A. Cucchieri, K. Langfeld, T. Mendes, Nucl. Phys. B687 (2004) 76–100.
  * Oliveira and Silva [2007] O. Oliveira, P. Silva, Braz.J. Phys. 37 (2007) 201–207.
  * Sternbeck et al. [2007] A. Sternbeck, K. Maltman, L. von Smekal, A. Williams, E. Ilgenfritz, M. Muller-Preussker, PoS LATTICE2007 (2007) 256\.
  * Boucaud et al. [2009] P. Boucaud, F. De Soto, J. Leroy, A. Le Yaouanc, J. Micheli, et al., Phys. Rev. D79 (2009) 014508.
  * Horak et al. [2023] J. Horak, J. M. Pawlowski, J. Turnwald, J. M. Urban, N. Wink, S. Zafeiropoulos, Phys. Rev. D 107 (2023) 076019\.
  * Cyrol et al. [2018] A. K. Cyrol, M. Mitter, J. M. Pawlowski, N. Strodthoff, Phys. Rev. D97 (2018) 054006\.
  * Corell et al. [2018] L. Corell, A. K. Cyrol, M. Mitter, J. M. Pawlowski, N. Strodthoff, SciPost Phys. 5 (2018) 066.
  * Deur et al. [2024] A. Deur, S. J. Brodsky, C. D. Roberts, Prog. Part. Nucl. Phys. 134 (2024) 104081.
  * de Teramond et al. [2024] G. F. de Teramond, A. Paul, S. J. Brodsky, A. Deur, H. G. Dosch, T. Liu, R. S. Sufian, arXiv:2403.16126 [hep-ph].
  * Cornwall [1982] J. M. Cornwall, Phys. Rev. D 26 (1982) 1453.
  * Aguilar and Papavassiliou [2010] A. C. Aguilar, J. Papavassiliou, Nucl. Phys. Proc. Suppl. 199 (2010) 172–177.
  * Aguilar et al. [2010] A. C. Aguilar, D. Binosi, J. Papavassiliou, J. High Energy Phys. 07 (2010) 002.
  * Boucaud et al. [2011] P. Boucaud, D. Dudal, J. Leroy, O. Pene, J. Rodriguez-Quintero, J. High Energy Phys. 12 (2011) 018.
  * von Smekal et al. [2009] L. von Smekal, K. Maltman, A. Sternbeck, Phys. Lett. B 681 (2009) 336–342.
  * Aguilar et al. [2009] A. C. Aguilar, D. Binosi, J. Papavassiliou, J. High Energy Phys. 11 (2009) 066.
  * Ilgenfritz et al. [2007] E.-M. Ilgenfritz, M. Muller-Preussker, A. Sternbeck, A. Schiller, I. Bogolubsky, Braz.J. Phys. 37 (2007) 193–200.
  * Bogolubsky et al. [2007] I. Bogolubsky, E. Ilgenfritz, M. Muller-Preussker, A. Sternbeck, PoS LATTICE2007 (2007) 290\.
  * Cucchieri and Mendes [2007] A. Cucchieri, T. Mendes, PoS LATTICE2007 (2007) 297\.
  * Oliveira and Silva [2012] O. Oliveira, P. J. Silva, Phys. Rev. D86 (2012) 114513\.
  * Braun et al. [2010] J. Braun, H. Gies, J. M. Pawlowski, Phys. Lett. B684 (2010) 262–267.
  * Fischer et al. [2009] C. S. Fischer, A. Maas, J. M. Pawlowski, Annals Phys. 324 (2009) 2408–2437.
  * Silva and Oliveira [2004] P. J. Silva, O. Oliveira, Nucl. Phys. (2004) 177.
  * Suman and Schilling [1996] H. Suman, K. Schilling, Phys. Lett. B 373 (1996) 314–318.
  * Cucchieri et al. [2018] A. Cucchieri, D. Dudal, T. Mendes, O. Oliveira, M. Roelfs, P. J. Silva, Phys. Rev. D 98 (2018) 091504\.
  * Davies et al. [1988] C. T. H. Davies, G. G. Batrouni, G. R. Katz, A. S. Kronfeld, G. P. Lepage, K. G. Wilson, P. Rossi, B. Svetitsky, Phys. Rev. D 37 (1988) 1581\.
  * Barrett et al. [1994] R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo, C. Romine, H. van der Vorst, Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods, Society for Industrial and Applied Mathematics, 1994.
  * Wilson [1974] K. G. Wilson, Phys. Rev. D10 (1974) 2445–2459.
  * Necco and Sommer [2002] S. Necco, R. Sommer, Nucl. Phys. B 622 (2002) 328–346.
  * Boucaud et al. [2018] P. Boucaud, F. De Soto, K. Raya, J. Rodriguez-Quintero, S. Zafeiropoulos, Phys. Rev. D98 (2018) 114515.
  * Duarte et al. [2016] A. G. Duarte, O. Oliveira, P. J. Silva, Phys. Rev. D 94 (2016) 014502.
  * Edwards and Joo [2005] R. G. Edwards, B. Joo (SciDAC, LHPC, UKQCD), Nucl. Phys. B Proc. Suppl. 140 (2005) 832\.
  * Pippig [2013] M. Pippig, SIAM Journal on Scientific Computing 35 (2013) C213–C236.
  * Oliveira et al. [2019] O. Oliveira, P. J. Silva, J.-I. Skullerud, A. Sternbeck, Phys. Rev. D 99 (2019) 094506\.
  * Kızılersü et al. [2021] A. Kızılersü, O. Oliveira, P. J. Silva, J.-I. Skullerud, A. Sternbeck, Phys. Rev. D 103 (2021) 114515.
  * Binosi et al. [2017] D. Binosi, C. D. Roberts, J. Rodriguez-Quintero, Phys. Rev. D 95 (2017) 114009.
  * Dudal et al. [2018] D. Dudal, O. Oliveira, P. J. Silva, Annals Phys. 397 (2018) 351–364.
  * Ferreira and Papavassiliou [2023] M. N. Ferreira, J. Papavassiliou, Particles 6 (2023) 312–363.
  * Aguilar et al. [2020] A. C. Aguilar, F. De Soto, M. N. Ferreira, J. Papavassiliou, J. Rodríguez-Quintero, S. Zafeiropoulos, Eur. Phys. J. C80 (2020) 154\.


# Determining Which Binary Component Hosts the TESS Transiting Planet

Kathryn V. Lester NASA Ames Research Center, Moffett Field, CA 94035, USA
Steve B. Howell NASA Ames Research Center, Moffett Field, CA 94035, USA  David
R. Ciardi NASA Exoplanet Science Institute, Caltech/IPAC, Pasadena, CA 91125,
USA  Rachel A. Matson U.S. Naval Observatory, Washington, D.C. 20392, USA
Kathryn Lester<EMAIL_ADDRESS>
(Accepted to AJ on )

###### Abstract

The NASA TESS mission has discovered many transiting planets orbiting bright
nearby stars, and high-resolution imaging studies have revealed that a number
of these exoplanet hosts reside in binary or multiple star systems. In such
systems, transit observations alone cannot determine which star in the binary
system actually hosts the orbiting planet. The knowledge of which star the
planet orbits is necessary for determining accurate physical properties for
the planet, especially its true radius and mean bulk density. We derived the
mean stellar densities for the components of 23 binary systems using the light
curve transit shape and the binary flux ratio from speckle imaging, then
tested the consistency with stellar models to determine which component is the
host star. We found that 70% of the TESS transiting planets in our sample
orbit the primary star.

††facilities: TESS, SOAR, Gemini North, Gemini South, WIYN

##  1 Introduction

The NASA TESS mission has discovered thousands of transiting planets orbiting
bright nearby stars (Ricker et al., 2014; Guerrero et al., 2021), and many
more planet candidates will be found as the TESS Extended Mission continues.
Follow-up high resolution imaging studies have revealed that a number of these
exoplanet host stars reside in binary or multiple star systems (e.g., Furlan
et al., 2017; Matson et al., 2018; Howell et al., 2021; Lester et al., 2021).
In these systems, transit observations alone cannot determine which binary
component star actually hosts the orbiting planet.

Identifying the correct host component is essential for accurate planet
properties. Correcting the transit depth for dilution from a companion greatly
affects the planet’s derived radius and thus the mean density and composition
(Ciardi et al., 2015; Furlan & Howell, 2017; Hirsch et al., 2017). Dilution
from a bright companion star also affects the host star properties determined
spectroscopically (Furlan & Howell, 2020), such as the stellar radius used to
calculate the planet radii. Lastly, binary components of different effective
temperatures have different habitable zones (HZ), so a transiting planet may
lie within the HZ if orbiting the primary star but outside the HZ if orbiting
the secondary star. Therefore, studies of exoplanet demographics and Earth-
like, habitable planets must take into account whether each planet hosts
resides in a multistar system (Savel et al., 2020), and if so, determine the
binary properties and which component actually hosts the planet.

Seager & Mallén-Ornelas (2003) derived a method to determine the stellar
density of a planet host star using the transit geometry, which can be
extended to planet hosts in multistar systems if the transit depth is
corrected for the flux contribution of each component. Payne et al. (2018)
used this method to identify the planet host star in 8 K2 binary systems and
found that the brighter primary star was more likely to host the planet in 5
out of 8 systems. They also found that systems with fainter companions were
more likely to have primary host stars.

In this work, we extend their study to 23 binary systems found to host
transiting planets by TESS. We introduce our sample in Section 2, then
describe the transit modeling and stellar density calculations in Section 3.
We identify the planet host components and present our results in Section 4,
and discuss the implications of our study in Section 5.

##  2 Sample

###  2.1 Binary TOI’s

We compiled a list of over 300 TOI exoplanet host candidates found to be
binary systems through speckle interferometry at Gemini North & South (Lester
et al., 2021), SOAR (Ziegler et al., 2020, 2021), and WIYN (Howell et al.,
2021). Systems suitable for this stellar density analysis must pass several
requirements by having:

  * •

binary separation $r<1.2$”, so the companions are well characterized by
speckle imaging

  * •

at least two consecutive transits, so the orbital period can be constrained
(Seager & Mallén-Ornelas, 2003)

  * •

flat bottom transits, so the contact times can be measured (Seager & Mallén-
Ornelas, 2003)

  * •

transits with sufficient signal-to-noise ratio (SNR $>3$) to model

  * •

planet candidate or confirmed planet status111as of December 2021 on
https://exofop.ipac.caltech.edu/tess/ on the Exoplanet Follow-up Observing
Program (ExoFOP) website to avoid background eclipsing binaries.

Of the 300+ TOI binaries evaluated, 106 systems had angular separations less
than 1.2”, and 23 systems passed all of the requirements for this analysis.
Our final sample is listed in Table 2.2 and includes 20 binary systems and 3
triple systems. Most of the rejected systems had V-shaped (grazing) transits
or were found to be false positives through follow-up photometry by the
community. Planets around M-type hosts stars may also create V-shaped transits
due to different limb darkening values than Solar type stars, but we had very
few M-type binaries in our sample to begin with.

###  2.2 Stellar Parameters

We adopted the effective temperatures ($T_{\rm eff}$) and uncertainties of the
primary components listed on ExoFOP from either the TESS Input Catalog (TIC,
Stassun et al., 2019) or a spectroscopic analysis (when available). We then
used magnitude difference ($\Delta m$) between the binary components from
speckle imaging to calculate the effective temperatures of the secondary and
tertiary components using a polynomial fit to the Modern Mean Dwarf Stellar
Color and Effective Temperature Sequence (Pecaut & Mamajek, 2013) as described
in Matson et al. (2019). For each primary star, we calculated the $I$-band
magnitude differences of all stars with cooler effective temperatures, fit a
seventh order polynomial to the these values as a function of $T_{\rm eff}$,
and interpolated the $T_{\rm eff}$ of the companions at the observed $\Delta
m$ values. We estimated the uncertainties in $T_{\rm eff}$ of the secondary or
tertiary stars based on the polynomial fit and the uncertainty in $\Delta m$,
then combined this value in quadrature with the uncertainties in the primary
star’s $T_{\rm eff}$. We adopted $\Delta m$ uncertainties of 0.2 mag (Howell
et al., 2021) for all systems, unless a specific value was listed on ExoFOP.
Metallicity measurements were not available for most systems, so we adopted
solar values where needed. The stellar properties of our binary TOI’s are
listed in Table 2.2, with the TIC number, the effective temperature of each
component, the metallicity, the binary angular separation(s), the magnitude
difference(s) from speckle imaging, and the corresponding references.

Table 1: Stellar Parameters

TOI | TIC |  $T_{\rm eff1}$ (K) |  $T_{\rm eff2}$ (K) |  $T_{\rm eff3}$ (K) | [Fe/H] |  $r_{12}$ (”) |  $r_{13}$ (”) |  $\Delta m_{12}$ (mag) |  $\Delta m_{13}$ (mag) | Ref.
---|---|---|---|---|---|---|---|---|---|---
141 | 403224672 | $5800\pm 110$ | $3640\pm 150$ | $3470\pm 150$ |  $\cdots$ | 0.49 | 1.33 | $4.63\pm 0.10$ | $5.60\pm 0.49$ | 1,2
159 | 394657039 | $6980\pm 150$ | $4760\pm 150$ |  $\cdots$ | $0.17\pm 0.12$ | 0.64 |  $\cdots$ | $3.30\pm 0.20$ |  $\cdots$ | 1,3
172 | 29857954 | $5760\pm 130$ | $3510\pm 150$ |  $\cdots$ | $0.15\pm 0.08$ | 1.13 |  $\cdots$ | $5.36\pm 0.20$ |  $\cdots$ | 2,4
194 | 211438925 | $5940\pm 150$ | $5320\pm 200$ |  $\cdots$ |  $\cdots$ | 0.27 |  $\cdots$ | $1.00\pm 0.20$ |  $\cdots$ | 1,3
264 | 122612091 | $6250\pm 100$ | $3859\pm 150$ |  $\cdots$ | $-0.06\pm 0.09$ | 0.65 |  $\cdots$ | $4.46\pm 0.20$ |  $\cdots$ | 2,5
295 | 327952677 | $5860\pm 120$ | $3800\pm 150$ |  $\cdots$ |  $\cdots$ | 0.98 |  $\cdots$ | $4.10\pm 0.20$ |  $\cdots$ | 1,3
309 | 228507250 | $5330\pm 250$ | $5140\pm 250$ |  $\cdots$ |  $\cdots$ | 0.34 |  $\cdots$ | $1.03\pm 0.20$ |  $\cdots$ | 1,2
322 | 198077394 | $6500\pm 200$ | $6080\pm 200$ |  $\cdots$ |  $\cdots$ | 0.18 |  $\cdots$ | $0.70\pm 0.20$ |  $\cdots$ | 1,3
462 | 420049884 | $5500\pm 150$ | $5240\pm 150$ |  $\cdots$ |  $\cdots$ | 0.17 |  $\cdots$ | $0.54\pm 0.20$ |  $\cdots$ | 3,7
487 | 31852980 | $5450\pm 130$ | $4685\pm 150$ |  $\cdots$ |  $\cdots$ | 0.52 |  $\cdots$ | $1.25\pm 0.20$ |  $\cdots$ | 1,2
489 | 455096220 | $6100\pm 100$ | $3730\pm 150$ |  $\cdots$ | $0.11\pm 0.08$ | 0.94 |  $\cdots$ | $4.70\pm 0.20$ |  $\cdots$ | 8,9,10
492 | 17746821 | $6280\pm 100$ | $4043\pm 150$ |  $\cdots$ | $0.18\pm 0.08$ | 0.67 |  $\cdots$ | $3.89\pm 0.20$ |  $\cdots$ | 2,9
601 | 141363913 | $6220\pm 140$ | $4370\pm 150$ |  $\cdots$ |  $\cdots$ | 0.15 |  $\cdots$ | $3.10\pm 0.20$ |  $\cdots$ | 1,3
640 | 147977348 | $6440\pm 130$ | $3880\pm 150$ |  $\cdots$ |  $\cdots$ | 0.25 |  $\cdots$ | $4.80\pm 0.20$ |  $\cdots$ | 1,3
697 | 77253676 | $5450\pm 130$ | $3468\pm 150$ | $3427\pm 150$ |  $\cdots$ | 1.17 | 1.27 | $4.95\pm 0.20$ | $5.09\pm 0.20$ | 1,2
931 | 206474443 | $5930\pm 150$ | $5640\pm 200$ |  $\cdots$ |  $\cdots$ | 0.12 |  $\cdots$ | $0.50\pm 0.10$ |  $\cdots$ | 1,2
1217 | 248092710 | $6310\pm 130$ | $3750\pm 150$ |  $\cdots$ |  $\cdots$ | 1.17 |  $\cdots$ | $5.02\pm 0.20$ |  $\cdots$ | 1,2
1267 | 158561566 | $6400\pm 100$ | $6100\pm 150$ |  $\cdots$ | $0.12\pm 0.06$ | 0.29 |  $\cdots$ | $0.59\pm 0.20$ |  $\cdots$ | 11,12
1287 | 352764091 | $5890\pm 130$ | $4025\pm 150$ |  $\cdots$ |  $\cdots$ | 0.10 |  $\cdots$ | $3.34\pm 0.20$ |  $\cdots$ | 1,12
1288 | 365733349 | $6180\pm 100$ | $4578\pm 150$ | $3523\pm 150$ |  $\cdots$ | 0.07 | 1.12 | $2.57\pm 0.20$ | $5.90\pm 0.20$ | 1,12
1307 | 236815160 | $5250\pm 200$ | $3800\pm 250$ |  $\cdots$ |  $\cdots$ | 0.33 |  $\cdots$ | $3.05\pm 0.20$ |  $\cdots$ | 7,12
1409 | 115771549 | $6110\pm 200$ | $4480\pm 200$ |  $\cdots$ |  $\cdots$ | 0.46 |  $\cdots$ | $2.63\pm 0.20$ |  $\cdots$ | 1,2
1740 | 174041208 | $5370\pm 150$ | $3710\pm 180$ |  $\cdots$ |  $\cdots$ | 0.24 |  $\cdots$ | $3.60\pm 0.20$ |  $\cdots$ | 1,2

References. — 1\. Stassun et al. (2019), 2\. Lester et al. (2021), 3\. Ziegler
et al. (2020), 4\. Rodriguez et al. (2019), 5. Gillon et al. (2013), 6.
Dumusque et al. (2019), 7. Spectroscopic follow-up analysis on ExoFOP, 8\.
Stassun et al. (2017), 9. Bonomo et al. (2017), 10\. Ziegler et al. (2021),
11\. Buchhave et al. (2011), 12. Howell et al. (2021)

##  3 Stellar Densities

###  3.1 TESS Photometry

We downloaded the 2 minute Pre-Search Data Conditioning Simple Aperture
Photometry (PDCSAP Jenkins et al., 2016) photometry for each TOI from the
Mikulski Archive for Space Telescopes (MAST), removed the data with poor
quality flags, then normalized each light curve by dividing by the mean count
value. Several systems (TOI 159, 309, 322, 487, 697, 1307, 1740) showed large
scale stellar variability that would hinder transit fitting. We adopted the
preliminary ephemerides from ExoFOP to mask out the planet transits for each
system, then used the celerite222https://github.com/exoplanet-dev/celerite2
Gaussian processes code (Foreman-Mackey et al., 2017; Foreman-Mackey, 2018) to
fit and remove the stellar variability. Figure 1 shows an example of the
celerite fit and flattened light curve for TOI 322. Finally, we binned each
light curve into 10 minute intervals and kept only the data points within
$10-12$ hours of mid-transit in order to save on computation time.

Figure 1: Removal of stellar variability signal for TOI 322. The top panel
shows the normalized TESS photometry from sector 32 (black points) and the
celerite fit to the out-of-transit variability (pink line). The bottom panel
shows the flattened, transit light curve on the same scale.

###  3.2 Transit Modeling

We fit the transit geometry using the EXOFAST (Eastman et al., 2013) online
tool333https://exoplanetarchive.ipac.caltech.edu/cgi-bin/ExoFAST/nph-exofast
from the NASA Exoplanet Science Institute (NExSci). We adopted the transit
solutions listed on ExoFOP as initial parameters and ran the MCMC modeling
mode to produce robust best-fit transit parameters and uncertainties. We held
the orbital eccentricity fixed to zero for all systems. The online tool does
not have a TESS bandpass, so we used the $I$-band option as it has the most
similar central wavelength. EXOFAST also incorporates limb darkening based on
the stellar effective temperatures from Section 2.2.

We did not include any radial velocities in the transit fitting procedure,
because our goal was to apply this method to as many TOI systems as possible
rather than completing an in-depth analysis on a smaller number of systems.
Therefore, our fits do not have the same precision as the combined
photometry/radial velocity fits from the literature, but are sufficient for
our stellar density calculations. Our best-fit transit geometry parameters are
listed in Table 2, including the planet period ($P$), transit depth
($\delta$), total transit duration ($t_{T}$), impact paramter ($b$), planet
radius ($R_{pl}$), and transit SNR. Because the host star identity is not yet
known, the transit depth and planet radius listed here have not been corrected
for the flux for the companion(s); we discuss the planet radius corrections in
Section 4. Figures 2 and 3 show the phase folded light curve and the best-fit
transit model for all of our systems.

Figure 2: Transit light curves for TOI 141 – 492. The phase folded TESS
photometry is shown as the blue circles and the best fit EXOFAST transit model
is shown as the solid black line.

Figure 3: Transit light curves for TOI 601 – 1740. The phase folded TESS photometry is shown as the blue circles and the best fit EXOFAST transit model is shown as the solid black line.  Table 2: Transit Parameters  Planet |  $P$ (d) |  $\delta$ (ppt)**Not corrected for flux dilution from the stellar companion(s). Radius corrections are discussed in Section 4. |  $t_{T}$ (d) | $b$ |  $R_{pl}$ ($R_{J}$)**Not corrected for flux dilution from the stellar companion(s). Radius corrections are discussed in Section 4. | SNR
---|---|---|---|---|---|---
TOI 141 b | $1.007933\pm 0.000173$ | $0.20\pm 0.01$ | $0.066\pm 0.003$ | $0.65\pm 0.04$ | $0.16\pm 0.01$ | 6.9
TOI 159 b | $3.762837\pm 0.000002$ | $10.00\pm 0.09$ | $0.169\pm 0.002$ | $0.28\pm 0.08$ | $1.80\pm 0.77$ | 87.4
TOI 172 b | $9.476893\pm 0.000022$ | $3.11\pm 0.07$ | $0.196\pm 0.002$ | $0.23\pm 0.10$ | $0.97\pm 0.03$ | 7.8
TOI 194 b | $4.899643\pm 0.000003$ | $9.76\pm 0.11$ | $0.140\pm 0.001$ | $0.70\pm 0.02$ | $1.29\pm 0.04$ | 63.9
TOI 264 b | $2.216735\pm 0.000001$ | $4.35\pm 0.06$ | $0.139\pm 0.003$ | $0.56\pm 0.01$ | $1.07\pm 0.05$ | 57.6
TOI 295 b | $3.397068\pm 0.000293$ | $10.20\pm 0.43$ | $0.099\pm 0.003$ | $0.80\pm 0.03$ | $1.36\pm 0.52$ | 18.8
TOI 309 b | $6.156681\pm 0.000807$ | $19.60\pm 5.71$ | $0.118\pm 0.006$ | $0.92\pm 0.32$ | $2.16\pm 3.44$ | 24.9
TOI 322 b | $3.990584\pm 0.000055$ | $3.33\pm 0.12$ | $0.154\pm 0.002$ | $0.75\pm 0.06$ | $1.03\pm 0.11$ | 27.3
TOI 462 b | $4.107456\pm 0.000027$ | $0.79\pm 0.10$ | $0.132\pm 0.008$ | $0.42\pm 0.33$ | $0.46\pm 0.24$ | 4.2
TOI 487 b | $7.414321\pm 0.000480$ | $0.35\pm 0.04$ | $0.158\pm 0.007$ | $0.39\pm 0.32$ | $0.33\pm 0.18$ | 2.9
TOI 489 b | $3.646659\pm 0.000004$ | $8.05\pm 0.27$ | $0.163\pm 0.003$ | $0.39\pm 0.22$ | $1.29\pm 0.15$ | 24.2
TOI 492 b | $3.122008\pm 0.000003$ | $6.19\pm 0.18$ | $0.153\pm 0.004$ | $0.73\pm 0.06$ | $1.31\pm 0.05$ | 27.1
TOI 601 b | $3.477689\pm 0.000708$ | $1.16\pm 0.19$ | $0.092\pm 0.008$ | $0.60\pm 0.41$ | $0.51\pm 0.25$ | 7.0
TOI 640 b | $5.003812\pm 0.000112$ | $7.70\pm 0.15$ | $0.152\pm 0.001$ | $0.89\pm 0.01$ | $1.78\pm 0.29$ | 55.4
TOI 697 b | $8.607870\pm 0.000055$ | $0.50\pm 0.08$ | $0.147\pm 0.007$ | $0.50\pm 0.34$ | $0.42\pm 0.23$ | 5.8
TOI 931 b | $3.245232\pm 0.000179$ | $6.21\pm 0.22$ | $0.116\pm 0.002$ | $0.28\pm 0.20$ | $0.78\pm 0.09$ | 18.8
TOI 1217 b | $41.485403\pm 0.000492$ | $0.95\pm 0.17$ | $0.256\pm 0.021$ | $0.68\pm 0.45$ | $0.89\pm 0.54$ | 7.1
TOI 1267 b | $6.790125\pm 0.000029$ | $2.10\pm 0.08$ | $0.247\pm 0.004$ | $0.31\pm 0.20$ | $0.89\pm 0.06$ | 7.1
TOI 1287 b | $9.597430\pm 0.000242$ | $0.42\pm 0.07$ | $0.214\pm 0.010$ | $0.61\pm 0.38$ | $0.34\pm 0.28$ | 4.7
TOI 1288 b | $2.699840\pm 0.000009$ | $2.73\pm 0.08$ | $0.109\pm 0.002$ | $0.82\pm 0.03$ | $0.86\pm 0.16$ | 41.9
TOI 1307 b | $2.534599\pm 0.000053$ | $6.19\pm 0.18$ | $0.100\pm 0.002$ | $0.87\pm 0.01$ | $1.16\pm 0.06$ | 37.1
TOI 1409 b | $4.717965\pm 0.000321$ | $3.00\pm 0.07$ | $0.192\pm 0.003$ | $0.87\pm 0.01$ | $1.51\pm 0.09$ | 33.3
TOI 1740 b | $19.428487\pm 0.000918$ | $1.01\pm 0.22$ | $0.154\pm 0.012$ | $0.58\pm 0.38$ | $0.26\pm 0.23$ | 3.1

###  3.3 Observed Stellar Densities

We calculated the mean stellar density for each component
($\rho_{1},\rho_{2},\rho_{3}$) following Equations 8-11 in Payne et al.
(2018), which correct the observed transit depth for the flux contamination
from the other components. For example, if the planet orbits the primary star
in a binary system, the primary’s stellar density is

$\rho_{1}=\frac{3\pi}{P^{2}G}\Bigg{(}\frac{\bigg{(}1+\sqrt{(1+\frac{f_{2}}{f_{1}})\delta}~{}\bigg{)}^{2}-b^{2}\big{[}1-\sin^{2}(t_{T}\pi/P)\big{]}}{\sin^{2}(t_{T}\pi/P)}\Bigg{)}^{3/2}$
(1)

where $P$, $\delta$, $t_{T}$, and $b$ are the transit parameters listed in the
previous section, $f_{2}/f_{1}$ is the binary flux ratio, and $G$ is the
gravitational constant. If the planet orbits the secondary star, the flux
ratio term is inverted and the secondary’s stellar density is

$\rho_{2}=\frac{3\pi}{P^{2}G}\Bigg{(}\frac{\bigg{(}1+\sqrt{(1+\frac{f_{1}}{f_{2}})\delta}~{}\bigg{)}^{2}-b^{2}\big{[}1-\sin^{2}(t_{T}\pi/P)\big{]}}{\sin^{2}(t_{T}\pi/P)}\Bigg{)}^{3/2}.$
(2)

The Gemini and WIYN speckle data were taken in the 832 nm filter and the SOAR
data were taken in the $I$-band filter. Both filters are similar to the TESS
bandpass so we calculated the flux ratio in the TESS bandpass directly from
the speckle magnitude difference.

To measure the stellar density uncertainties, we varied all input parameters
within their uncertainties over 10,000 iterations, each time calculating the
stellar density if the planet orbits the primary, secondary, and tertiary (if
applicable). We then created a histogram of the resulting densities for each
component as shown in Figure 4, fit a Gaussian to each histogram, and adopted
the width as the $1\sigma$ uncertainties for each component. Our results are
listed in Table 3 with the stellar density and uncertainty of each binary
component. We found the biggest source of uncertainty in the stellar densities
to be the transit impact parameter, especially if the transit is near grazing
and the impact parameter’s uncertainties extend beyond $b=1$.

Figure 4: Example density histograms for TOI 172 from 10,000 iterations. The
histograms for the primary and secondary component are shown in blue and
green, respectively. The Gaussian fits, used to determine the $1\sigma$
uncertainties, are shown in solid black.

##  4 Results

###  4.1 Determining Host Component

We compared the observed stellar densities to model values to determine which
component provides the best match. To build our model density-$T_{\rm eff}$
curve, we started with TOI’s found to be single stars with Gemini speckle
interferometry (Lester et al., 2021). Speckle with Gemini is more sensitive to
companions than with WIYN or SOAR due to the larger aperture and better
angular resolution, so this sample would be the least contaminated by
unresolved binaries. We pulled the TIC parameters for each TOI, kept only
stars with uncertainties in radius, temperature, and mass less than 20%, and
calculated the mean density for each system. Following the method of Payne et
al. (2018), we then created a model density-$T_{\rm eff}$ curve by fitting a
6th-order polynomial to these temperatures and densities. Figures 5 and 6 show
the stellar density vs. effective temperature plots for all of the TOI’s in
our sample.

Next, we calculated the distance of the observed stellar parameters from this
density-$T_{\rm eff}$ model (in terms of $\sigma$) for each of the 10,000
iterations used in the previous section. Finally, we tallied how many times
each component was closest to the model and adopted this percentage as the
host star probability ($P_{1,2,3}$). We classified each system system as
“likely primary host” for $P_{1}\geq 0.7$, “likely secondary host” for
$P_{2}\geq 0.7$, “likely tertiary host” for $P_{3}\geq 0.7$, and “uncertain
host” for $0.3<P_{1}<0.7$. The host probability of each binary component and
classifications are listed in Table 3, and we discuss each system in detail in
the next section.

We also calculated the planet radius correction factor ($X_{R}$) using the
equations in Ciardi et al. (2015) for the systems where the host component was
identified (listed in Table 3). The radius correction for planets orbiting the
primary star depends only on the binary magnitude difference, while the radius
corrections for planets orbiting the secondary star depends on the magnitude
difference and the radius ratio of the stars. We do not have individual colors
or age information for our binaries, so we estimated the radii of each
component by interpolating the values in the Modern Mean Dwarf Stellar Color
and Effective Temperature Sequence (Pecaut & Mamajek, 2013) at the effective
temperatures of each component. As a self check, we also estimated the radius
ratio using a 1 Gyr MIST isochrone (Dotter, 2016; Choi et al., 2016) and found
similar results to within 5%. The planet radius corrections are typically
small ($\overline{X_{R}}=1.07$) when the planet orbits the primary star but
large ($\overline{X_{R}}=4.97$) when the planet orbits the secondary star.

Figure 5: Stellar density vs. effective temperature plots for TOI 141 – 492.
The dervied value for the primary component is shown as the blue circle, the
secondary component is shown as the green diamond, and any tertiary components
are shown as red triangles. The solid purple line shows the model
density-$T_{\rm eff}$ relation built from the single TOI’s (purple points).
The host probabilities of each component are also listed.

Figure 6: Stellar density vs. effective temperature plots for TOI 601 – 1740. The dervied value for the primary component is shown as the blue circle, the secondary component is shown as the green diamond, and any tertiary components are shown as red triangles. The solid purple line shows the model density-$T_{\rm eff}$ relation built from the single TOI’s (purple points). The host probabilities of each component are also listed.  Table 3: Stellar Densities & Host Probabilities  TOI |  $\rho_{1}$ (cgs) |  $\rho_{2}$ (cgs) |  $\rho_{3}$ (cgs) |  $P_{1}$ |  $P_{2}$ |  $P_{3}$ | Result | $X_{R}$
---|---|---|---|---|---|---|---|---
141 | $1.07\pm 0.19$ | $1.70\pm 0.27$ | $2.18\pm 0.46$ | 0.95 | 0.03 | 0.02 | Likely primary host | 1.01
159 | $0.58\pm 0.04$ | $1.43\pm 0.15$ |  $\cdots$ | 0.42 | 0.58 |  $\cdots$ | Uncertain host |  $\cdots$
172 | $0.84\pm 0.06$ | $3.38\pm 0.40$ |  $\cdots$ | 0.96 | 0.04 |  $\cdots$ | Likely primary host | 1.01
194 | $0.73\pm 0.04$ | $0.97\pm 0.06$ |  $\cdots$ | 0.72 | 0.28 |  $\cdots$ | Likely primary host | 1.18
264 | $0.40\pm 0.03$ | $1.46\pm 0.19$ |  $\cdots$ | 1.00 | 0.00 |  $\cdots$ | Likely primary host | 1.01
295 | $0.96\pm 0.16$ | $6.83\pm 1.21$ |  $\cdots$ | 1.00 | 0.00 |  $\cdots$ | Likely primary host | 1.01
309 | $2.38\pm 3.47$ | $3.62\pm 7.09$ |  $\cdots$ | 0.31 | 0.69 |  $\cdots$ | Uncertain host |  $\cdots$
322 | $0.31\pm 0.07$ | $0.36\pm 0.08$ |  $\cdots$ | 0.80 | 0.20 |  $\cdots$ | Likely primary host | 1.23
462 | $1.06\pm 0.32$ | $1.09\pm 0.32$ |  $\cdots$ | 0.73 | 0.27 |  $\cdots$ | Likely primary host | 1.27
487 | $1.08\pm 0.29$ | $1.16\pm 0.29$ |  $\cdots$ | 0.93 | 0.07 |  $\cdots$ | Likely primary host | 1.15
489 | $0.62\pm 0.21$ | $2.60\pm 0.57$ |  $\cdots$ | 0.88 | 0.12 |  $\cdots$ | Likely primary host | 1.01
492 | $0.28\pm 0.06$ | $1.16\pm 0.18$ |  $\cdots$ | 1.00 | 0.00 |  $\cdots$ | Likely primary host | 1.01
601 | $1.98\pm 1.53$ | $2.56\pm 0.91$ |  $\cdots$ | 0.19 | 0.81 |  $\cdots$ | Likely secondary host | 8.08
640 | $0.22\pm 0.02$ | $3.37\pm 0.56$ |  $\cdots$ | 0.90 | 0.10 |  $\cdots$ | Likely primary host | 1.01
697 | $1.45\pm 0.62$ | $2.41\pm 0.67$ | $2.52\pm 0.65$ | 0.87 | 0.10 | 0.02 | Likely primary host | 1.01
931 | $1.60\pm 0.33$ | $1.73\pm 0.35$ |  $\cdots$ | 0.19 | 0.81 |  $\cdots$ | Likely secondary host | 1.86
1217 | $0.92\pm 1.30$ | $2.13\pm 0.82$ |  $\cdots$ | 0.46 | 0.54 |  $\cdots$ | Uncertain host |  $\cdots$
1267 | $0.32\pm 0.08$ | $0.32\pm 0.08$ |  $\cdots$ | 0.80 | 0.20 |  $\cdots$ | Likely primary host | 1.31
1287 | $0.52\pm 0.34$ | $0.69\pm 0.30$ |  $\cdots$ | 0.84 | 0.16 |  $\cdots$ | Likely primary host | 1.02
1288 | $0.39\pm 0.07$ | $0.81\pm 0.11$ | $5.69\pm 0.94$ | 0.95 | 0.00 | 0.05 | Likely primary host | 1.05
1307 | $0.42\pm 0.04$ | $1.58\pm 0.20$ |  $\cdots$ | 0.62 | 0.38 |  $\cdots$ | Uncertain host |  $\cdots$
1409 | $0.10\pm 0.01$ | $0.23\pm 0.03$ |  $\cdots$ | 1.00 | 0.00 |  $\cdots$ | Likely primary host | 1.04
1740 | $2.00\pm 1.56$ | $3.47\pm 1.78$ |  $\cdots$ | 0.41 | 0.59 |  $\cdots$ | Uncertain host |  $\cdots$

###  4.2 Notes on Individual Systems

TOI 141 – The observed parameters of the primary component match the model
relation well; TOI 141 has a 95% probability that the primary star hosts the
transiting planet. This is consistent with past studies of TOI 141, as
Espinoza et al. (2020) confirmed the planetary nature of TOI 141 b and found
that it orbits the primary star through spectroscopic follow-up. Because TOI
141 b likely orbits the primary star, the radius correction factor is only
1.01.

TOI 159 – The observed parameters of both components match the model relation
equally well, so TOI 159 has a 42% and 58% probability that the primary and
secondary star host the transiting planet, respectively. The transit geometry
and stellar densities are well constrained, so perhaps the binary magnitude
difference is causing a poor fit to the model relation. TOI 159 b creates a
deep transit, so it’s also possible that it could be a red dwarf companion
instead. This would be a good candidate to follow-up to confirm the binary
parameters and nature of the transiting object. This would also help constrain
the planet candidate’s radius, as the correction factor for TOI 159 b would be
1.02 if it orbits the primary and 10.66 if it orbits the secondary. However,
such a large correction would cause the planet density to be unrealistically
small (Furlan & Howell, 2020).

TOI 172 – The observed parameters of the primary component match the model
relation well; TOI 172 has a 96% probability that the primary star hosts the
transiting planet. This is consistent with the past study confirming TOI 172 b
as a planet by Rodriguez et al. (2019) through spectroscopic follow-up.
Because TOI 141 b likely orbits the primary star, the radius correction factor
is only 1.01.

TOI 194 – TOI 194 has a 72% and 28% probability that the primary and secondary
star host the transiting planet, respectively. A primary component host star
is consistent with the study confirming TOI 194 b as a planet (WASP-20b) by
Anderson et al. (2015) through radial velocity analysis, but this study did
not list a stellar companion so the derived planet radius is underestimated by
a factor of 1.18.

TOI 264 – The observed parameters of the primary component match the model
relation well; TOI 264 has a 100% probability that the primary star hosts the
transiting planet. This is consistent with the past study confirming TOI 264 b
as a planet (WASP-72b) by Gillon et al. (2013). This study did not list a
stellar companion, but the radius correction factor is only 1.01, so the
planet radius would not change within the uncertainties.

TOI 295 – The observed parameters of the primary component match the model
relation well; TOI 295 has a 100% probability that the primary star hosts the
transiting planet. This is consistent with the past spectroscopic analysis to
confirm TOI 295 b as a planet (WASP-146b) by Anderson et al. (2018). This
study did not list a stellar companion, but the radius correction factor is
only 1.01, so the planet radius would not change within the uncertainties.

TOI 309 – Due to the large impact parameter found by the transit fit, the
stellar densities are poorly constrained. We found that TOI 309 has a 31% and
69% probability that the primary and secondary star host the transiting
planet, respectively, so the host star is uncertain. If the planet orbits the
primary, the radius correction factor would be 1.18, but if the planet orbits
the secondary, the correction factor would be 2.20 and largely change the
derived planet density and classification. Additional follow-up is needed to
rule out any background eclipsing binaries or false positive scenarios and to
better characterize the system before a host star can be identified.

TOI 322 – The observed parameters of the primary component match the model
relation slightly better than that of the secondary component. We found that
TOI 322 has an 80% and 20% probability that the primary and secondary star
host the transiting planet, respectively, so the radius correction would be a
factor of 1.23.

TOI 462 – The observed parameters of the primary component match the model
relation slightly better than that of the secondary component. We found that
TOI 462 has an 73% and 27% probability that the primary and secondary star
host the transiting planet, respectively, so the planet radius correction
would be a factor of 1.27.

TOI 487 – The observed parameters of the primary component match the model
relation well; TOI 487 has a 93% and 7% probability that the primary and
secondary star host the transiting planet, respectively. The radius correction
factor for TOI 487 b would be 1.15.

TOI 489 – The observed parameters of the primary component match the model
relation well; TOI 489 has an 88% and 12% probability that the primary and
secondary star host the transiting planet, respectively. This is consistent
with the past study confirming TOI 489 b as a planet (HAT-P-35b) by Bakos et
al. (2012) through spectroscopic analysis. This study did not list a stellar
companion, but the radius correction factor is only 1.01, so the planet radius
would not change within the uncertainties.

TOI 492 – The observed parameters of the primary component match the model
relation well; TOI 492 has a 100% probability that the primary star hosts the
transiting planet. This is consistent with the past study confirming TOI 492 b
as a planet (HAT-P-50b) by Hartman et al. (2015). This study did not list a
stellar companion, but did test light curve fits that included an unseen
companion and could not rule out scenarios with a faint secondary star (as
found later by speckle). However, the planet radius correction factor is only
1.01, so the derived planet radius would not change within the uncertainties.

TOI 601 – The observed parameters of the secondary component match the model
relation better than those of the primary component, making TOI 601 one of the
few “likely secondary host” systems in our sample. TOI 601 has a 19% and 81%
probability that the primary and secondary stars host the transiting planet,
respectively. However, the stellar density uncertainties are rather large due
to the low SNR transit, so we recommend additional follow-up to better
characterize this system, constrain the stellar densities, and confirm the
identity of the host star. This would also help constrain the planet radius,
as the correction factor for TOI 601 b would be minor (1.03) if it orbits the
primary but quite large (8.08) if it orbits the secondary. However, such a
large correction would cause the planet density to be unphysically small
(Furlan & Howell, 2020).

TOI 640 – The observed parameters of the primary component match the model
relation well; TOI 640 has a 90% probability that the primary star hosts the
transiting planet. TOI 640 b was confirmed to be a planet by Rodriguez et al.
(2021), but they reported that no stellar companions were found by SOAR, even
though the source paper did report a 0.25” companion (Ziegler et al., 2020).
Nonetheless, the secondary star is 4.8 mag fainter than the primary and
unlikely to be seen in spectroscopic follow-up, so the radial velocity signal
detected by Rodriguez et al. (2021) would also confirm the primary star hosts
the planet. The radius correction factor for TOI 640 b is only 1.01, so the
planet radius would not be largely affected.

TOI 697 – The observed parameters of the primary component match the model
relation well. TOI 697 has a 87%, 10%, and 2% probability that the primary,
secondary, and tertiary stars host the transiting planet, respectively.
Because the planet likely orbits the primary, the radius correction is a
factor of only 1.01.

TOI 931 – Neither the primary or secondary components fit the model relation
very well, but the secondary matches slightly better. We found that TOI 931
has a 19% and 81% probability that the primary and secondary stars host the
transiting planet, respectively. The radius correction factor for TOI 931 b
would be 1.86.

TOI 1217 – The observed parameters of both components match the model relation
due to the large uncertainties. We found that TOI 1217 has a 46% and 54%
probability that the primary and secondary stars host the transiting planet,
respectively, so the host star is uncertain until additional photometric
observations can better constrain the transit geometry. This would also help
constrain the planet radius, as the correction factor for TOI 1217 b would be
minor (1.01) if it orbits the primary but very large (26.99). However, such a
large correction would cause the planet density to be unphysically small
(Furlan & Howell, 2020).

TOI 1267 – The observed parameters of the primary component match the model
relation slightly better than those of the secondary. We found that TOI 1267
has a 80% and 20% probability that the primary and secondary stars host the
transiting planet, respectively. This is consistent with the study confirming
TOI 492 b as a planet (Kepler-14b) through spectroscopic follow-up by Buchhave
et al. (2011), who concluded the planet orbits the primary star using a
centroid analysis of the Kepler data and corrected the derived planet
parameters accordingly.

TOI 1287 – The observed parameters of the primary component match the model
relation well; TOI 1287 has an 84% and 16% probability that the primary and
secondary stars host the transiting planet, respectively. The radius
correction for TOI 1287 b would therefore be a factor of 1.02.

TOI 1288 – The observed parameters of the primary component match the model
relation well. TOI 1288 has a 95%, 0%, and 5% probability that the primary,
secondary, and tertiary stars host the transiting planet, respectively. The
radius correction for TOI 1288 b would therefore be a factor of 1.05.

TOI 1307 – Neither the primary or secondary components fit the model relation
very well, so the host star is uncertain. We found that TOI 1307 has a 62% and
38% probability that the primary and secondary stars host the transiting
planet, respectively. If the planet orbits the primary star, the radius
correction factor would be only 1.03, but if it orbits the secondary star, the
correction would be a factor of 6.54. The transit geometry is well
constrained, so the poor fit to the model relation is likely due to the
primary star’s evolutionary state. The ExoFOP website notes that TOI 1307 is
slightly evolved, which would cause the observed density to be lower than the
model (made from mostly main sequence TOI’s).

TOI 1409 – The observed parameters of the primary component match the model
relation well; we found that TOI 1409 has a 100% probability that the primary
stars hosts the transiting planet. The planet radius correction would
therefore be a factor of 1.04.

TOI 1740 – The observed parameters of both components match the model relation
due to the large uncertainties in the density, so the host star is uncertain.
We found that TOI 1740 has a 41% and 59% probability that the primary and
secondary stars host the transiting planet, respectively. If the planet orbits
the primary star, the radius correction factor would be only 1.02, but if it
orbits the secondary star, the correction would be a factor of 9.37.
Additional photometric observations are needed to better constrain the transit
geometry and identify the planet host.

###  4.3 Overall Results

Figure 7: The probability of a primary star host ($P_{1}$) vs. the binary
magnitude difference ($\Delta m$) for each binary in our sample. The “likely
primary host” systems are plotted in blue, the “likely secondary host” systems
are plotted in green, and the “uncertain host” systems are plotted in grey.
The preference for primary star hosts in systems with faint companions is
easily visible.

Here are the results of our study:

  * •

Likely primary host – 16 out of 23 systems (70%)

  * •

Likely secondary host – 2 out of 23 systems (9%)

  * •

Likely tertiary host – 0 out of 3 systems (0%)

  * •

Uncertain host – 5 out of 23 systems (21%)

We found that most transiting planets in our sample (70%) orbit the primary
star in their binary system. Figure 7 shows the primary host star probability
as a function of the the binary magnitude difference. While we do not see the
strong trend found by Payne et al. (2018), where the host component in equal
brightness pairs were uncertain due to the similarity in temperature and
density, the equal brightness systems in our sample do have slightly lower
primary host probabilities that unequal brightness pairs. Furthermore, the
systems with faint secondaries ($\Delta m\geq$ 4 mag) are even more likely to
have a primary star host (89%), as seen in Figure 7. This is because the light
from the primary star washes out the transit signal from a planet orbiting the
secondary (or tertiary) star until it becomes too difficult to detect (e.g.,
Lester et al., 2021). This result is consistent with that of Furlan & Howell
(2017), who found that faint secondary host stars would result in planet
densities that are too low to be physically realistic. Finally, we did not
find any trends with planet radius or binary separation due to the small
number of secondary host systems.

Note that our results apply to transiting planets in binary systems, not to
all planets. The main selection effects of our study are that speckle imaging
cannot detect very close binary companions ($<$20 mas) and that our transit
signal-to-noise limit may exclude small planets. However, past studies (e.g.,
Kraus et al., 2016; Moe & Kratter, 2021; Howell et al., 2021) have found that
binary host separation peaks around 100 AU, which is well within the speckle
angular resolution limits for most TESS stars, and that the number of
companions rapidly decreases towards small separations. Next, the transits of
small planets ($\leq 2R_{E}$) in binary systems are not seen by TESS due to
the observational bias from third light contamination (Lester et al., 2021),
but the transits of giant planets commonly detected with TESS can still be
seen in multistar systems. While a full completeness analysis is beyond the
scope of this paper, our observational sample broadly represents the binary
planet hosts found by TESS.

##  5 Conclusions

Using the transit light curve from TESS and binary parameters from speckle
imaging, we determined the stellar densities of each component in 23 binary
exoplanet host systems. By comparing these observed densities with models, we
determined which binary component was more likely to host the planet. We found
that 70% of the TESS transiting planets in our sample are likely to orbit the
primary star in their binary host system. If this holds for all transit
observations, it would be reasonable in future studies to assume that the
primary star hosts the transiting planet when no additional observations are
available. The planets in these systems would still need a small radius
correction, as the fainter companions do slightly dilute the planet transit
(Ciardi et al., 2015).

However, we found several systems have a non-negligible chance that the planet
orbits the secondary star, which would greatly affect the resulting planet
parameters. Further follow-up is recommended for systems with an uncertain
host star or poor density fits – likely due to low SNR transits, small number
of transits, or small binary $\Delta m$. This could include additional transit
observations with TESS or groundbased photometry, spectroscopic observations
to search for the planet’s radial velocity signal, or multi-wavelength
photometric follow-up to further refine the stellar effective temperatures.

Identifying which binary component hosts the transiting planet is essential
for determining the correct transit depth, the true planet radius and mean
density, and the host star’s properties and habitable zone (Ciardi et al.,
2015; Hirsch et al., 2017; Furlan & Howell, 2017, 2020). Studies of exoplanet
demographics must take into account whether each planet host resides in a
multi-star system to ensure the derived planet properties and occurrence rates
are corrected accordingly (Bouma et al., 2018; Savel et al., 2020). Knowing
which planet hosts are binary systems and which component hosts the planet is
also important for target selection of upcoming coronagraph observations
(e.g., by Roman Space Telescope or LUVOIR). In order to directly image planets
in well-resolved binary systems, one must know ahead of time which star to
block with the coronagraph or if it is possible to block the light from both
stars (Bendek et al., 2021).

We would like to thank the anonymous referee for their helpful comments. KVL
is supported by an appointment to the NASA Postdoctoral Program at the NASA
Ames Research Center, administered by Oak Ridge Associated Universities under
contract with NASA. This paper includes data collected with the TESS mission
obtained from the MAST data archive at the Space Telescope Science Institute
(STScI). Funding for the TESS mission is provided by the NASA Explorer
Program. STScI is operated by the Association of Universities for Research in
Astronomy, Inc., under NASA contract NAS 5–26555. This work also made use of
the Exoplanet Follow-up Observation Program website (ExoFOP, 2019) and EXOFAST
(Eastman et al., 2013) tool provided by the NASA Exoplanet Archive, which is
operated by the California Institute of Technology under contract with the
National Aeronautics and Space Administration under the Exoplanet Exploration
Program.

## References

  * (1)
  * Anderson et al. (2015) Anderson, D. R., Collier Cameron, A., Hellier, C., et al. 2015, A&A, 575, A61. doi:10.1051/0004-6361/201423591
  * Anderson et al. (2018) Anderson, D. R., Bouchy, F., Brown, D. J. A., et al. 2018, arXiv:1812.09264
  * Astropy Collaboration et al. (2013) Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33. doi:10.1051/0004-6361/201322068
  * Astropy Collaboration et al. (2018) Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123. doi:10.3847/1538-3881/aabc4f
  * Bakos et al. (2012) Bakos, G. Á., Hartman, J. D., Torres, G., et al. 2012, AJ, 144, 19. doi:10.1088/0004-6256/144/1/19
  * Bendek et al. (2021) Bendek, E. A., Belikov, R., Sirbu, D., et al. 2021, Proc. SPIE, 11823, 1182311. doi:10.1117/12.2594992
  * Bonomo et al. (2017) Bonomo, A. S., Desidera, S., Benatti, S., et al. 2017, A&A, 602, A107. doi:10.1051/0004-6361/201629882
  * Bouma et al. (2018) Bouma, L. G., Masuda, K., & Winn, J. N. 2018, AJ, 155, 244. doi:10.3847/1538-3881/aabfb8
  * Buchhave et al. (2011) Buchhave, L. A., Latham, D. W., Carter, J. A., et al. 2011, ApJS, 197, 3. doi:10.1088/0067-0049/197/1/3
  * Choi et al. (2016) Choi, J., Dotter, A., Conroy, C., et al. 2016, ApJ, 823, 102. doi:10.3847/0004-637X/823/2/102
  * Ciardi et al. (2015) Ciardi, D. R., Beichman, C. A., Horch, E. P., et al. 2015, ApJ, 805, 16
  * Dotter (2016) Dotter, A. 2016, ApJS, 222, 8. doi:10.3847/0067-0049/222/1/8
  * Dumusque et al. (2019) Dumusque, X., Turner, O., Dorn, C., et al. 2019, A&A, 627, A43. doi:10.1051/0004-6361/201935457
  * Eastman et al. (2013) Eastman, J., Gaudi, B. S., & Agol, E. 2013, PASP, 125, 83. doi:10.1086/669497
  * Espinoza et al. (2020) Espinoza, N., Brahm, R., Henning, T., et al. 2020, MNRAS, 491, 2982. doi:10.1093/mnras/stz3150
  * ExoFOP (2019) ExoFOP. 2019. Exoplanet Follow-up Observing Program - TESS. IPAC. doi: 10.26134/EXOFOP3
  * Foreman-Mackey et al. (2017) Foreman-Mackey, D., Agol, E., Ambikasaran, S., et al. 2017, AJ, 154, 220. doi:10.3847/1538-3881/aa9332
  * Foreman-Mackey (2018) Foreman-Mackey, D. 2018, Research Notes of the American Astronomical Society, 2, 31. doi:10.3847/2515-5172/aaaf6c
  * Furlan et al. (2017) Furlan, E., Ciardi, D. R., Everett, M. E., et al. 2017, AJ, 153, 71. doi:10.3847/1538-3881/153/2/71
  * Furlan & Howell (2017) Furlan, E. & Howell, S. B. 2017, AJ, 154, 66. doi:10.3847/1538-3881/aa7b70
  * Furlan & Howell (2020) Furlan, E. & Howell, S. B. 2020, ApJ, 898, 47. doi:10.3847/1538-4357/ab9c9c
  * Gillon et al. (2013) Gillon, M., Anderson, D. R., Collier-Cameron, A., et al. 2013, A&A, 552, A82. doi:10.1051/0004-6361/201220561
  * Guerrero et al. (2021) Guerrero, N. M., Seager, S., Huang, C. X., et al. 2021, ApJS, 254, 39. doi:10.3847/1538-4365/abefe1
  * Harris et al. (2020) Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Nature, 585, 357. doi:10.1038/s41586-020-2649-2
  * Hartman et al. (2015) Hartman, J. D., Bhatti, W., Bakos, G. Á., et al. 2015, AJ, 150, 168. doi:10.1088/0004-6256/150/6/168
  * Hirsch et al. (2017) Hirsch, L. A., Ciardi, D. R., Howard, A. W., et al. 2017, AJ, 153, 117. doi:10.3847/1538-3881/153/3/117
  * Howell et al. (2021) Howell, S. B., Matson, R. A., Ciardi, D. R., et al. 2021, AJ, 161, 164. doi:10.3847/1538-3881/abdec6
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90. doi:10.1109/MCSE.2007.55
  * Jenkins et al. (2016) Jenkins, J. M., Twicken, J. D., McCauliff, S., et al. 2016, Proc. SPIE, 9913, 99133E. doi:10.1117/12.2233418
  * Kraus et al. (2016) Kraus, A. L., Ireland, M. J., Huber, D., et al. 2016, AJ, 152, 8
  * Lester et al. (2021) Lester, K. V., Matson, R. A., Howell, S. B., et al. 2021, AJ, 162, 75. doi:10.3847/1538-3881/ac0d06
  * Matson et al. (2018) Matson, R. A., Howell, S. B., Horch, E. P., et al. 2018, AJ, 156, 31
  * Matson et al. (2019) Matson, R. A., Howell, S. B., & Ciardi, D. R. 2019, AJ, 157, 211
  * Moe & Kratter (2021) Moe, M. & Kratter, K. M. 2021, MNRAS, 507, 3593. doi:10.1093/mnras/stab2328
  * Payne et al. (2018) Payne, A. N., Ciardi, D. R., Kane, S. R., et al. 2018, AJ, 156, 209. doi:10.3847/1538-3881/aae310
  * Pecaut & Mamajek (2013) Pecaut, M. J. & Mamajek, E. E. 2013, ApJS, 208, 9
  * Ricker et al. (2014) Ricker, G. R., Winn, J. N., Vanderspek, R., et al. 2014, Proc. SPIE, 9143, 914320. doi:10.1117/12.2063489
  * Rodriguez et al. (2019) Rodriguez, J. E., Quinn, S. N., Huang, C. X., et al. 2019, AJ, 157, 191. doi:10.3847/1538-3881/ab11d9
  * Rodriguez et al. (2021) Rodriguez, J. E., Quinn, S. N., Zhou, G., et al. 2021, AJ, 161, 194. doi:10.3847/1538-3881/abe38a
  * Savel et al. (2020) Savel, A. B., Dressing, C. D., Hirsch, L. A., et al. 2020, AJ, 160, 287. doi:10.3847/1538-3881/abc47d
  * Seager & Mallén-Ornelas (2003) Seager, S. & Mallén-Ornelas, G. 2003, ApJ, 585, 1038. doi:10.1086/346105
  * Stassun et al. (2017) Stassun, K. G., Collins, K. A., & Gaudi, B. S. 2017, AJ, 153, 136. doi:10.3847/1538-3881/aa5df3
  * Stassun et al. (2019) Stassun, K. G., Oelkers, R. J., Paegert, M., et al. 2019, AJ, 158, 138. doi:10.3847/1538-3881/ab3467
  * Virtanen et al. (2020) Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nature Methods, 17, 261. doi:10.1038/s41592-019-0686-2
  * Ziegler et al. (2020) Ziegler, C., Tokovinin, A., Briceño, C., et al. 2020, AJ, 159, 19
  * Ziegler et al. (2021) Ziegler, C., Tokovinin, A., Latiolais, M., et al. 2021, arXiv:2103.12076