aainstitutetext: Asia Pacific Center for Theoretical Physics, Postech, Pohang
37673, Koreabbinstitutetext: Department of Physics, Yonsei University, Seoul
03722, Koreaccinstitutetext: Department of Physics, Postech, Pohang 37673,
Korea

# Binary Black Holes and
Quantum Off-Shell Recursion

Kyoungho Cho b Kwangeon Kim a,c Kanghoon Lee<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The quantum off-shell recursion provides an efficient and universal
computational tool for loop-level scattering amplitudes. In this work, we
present a new comprehensive computational framework based on the quantum off-
shell recursion for binary black hole systems. Using the quantum perturbiner
method, we derive the recursions and solve them explicitly up to two-loop
order. We develop a power-counting prescription that enables the
straightforward separation of classical diagrams. We also devise a
classification scheme that optimizes the integration by parts (IBP) reduction
process, which makes higher-loop calculations more tractable. By employing the
soft expansion technique, we remove irrelevant terms from the loop integrands
and express them in terms of master integrals. We classify the one-loop and
the two-loop classical diagrams, and their loop integrands are represented by
linear combinations of the master integrals. Finally, we explicitly calculate
the classical scalar $2\to 2$ amplitudes in the potential region up to the 3PM
order and reproduce the known results.

††preprint: APCTP Pre2023 - 012

##  1 Introduction

Recently, quantum field theory (QFT) techniques have been employed to develop
a rigorous theoretical framework for studying binary black hole systems
despite the significant difference in the scale between elementary particles
and black holes. There has been substantial progress in this direction and
established a comprehensive theoretical framework that can accurately describe
the dynamics of binary black hole systems, especially in the inspiral phase
where the black holes are sufficiently separated. In this limit, black holes
can be considered as point particles up to a particular scale, allowing their
interactions to be treated as perturbations. In this context, the field theory
method Neill:2013wsa ; Luna:2017dtq ; Kosower:2018adc ; Foffa:2016rgu ;
Bjerrum-Bohr:2018xdl ; Chung:2018kqs ; Cheung:2018wkq ; Plefka:2018dpa ;
Bern:2019nnu ; Bern:2019crd ; Blumlein:2019zku ; Brandhuber:2019qpg ;
Cristofoli:2019ewu ; Cristofoli:2019neg ; Damour:2019lcq ; Foffa:2019hrb ;
Kalin:2019rwq ; KoemansCollado:2019ggb ; Bern:2020gjj ; Blumlein:2020znm ;
Bini:2020wpo ; Bini:2020uiq ; Bini:2020rzn ; Cheung:2020gyp ; Cheung:2020sdj ;
Cristofoli:2020uzm ; DiVecchia:2020ymx ; Kalin:2020mvi ; Kalin:2020fhe ;
Parra-Martinez:2020dzs ; Bern:2021yeh ; Bern:2021dqo ; Bjerrum-Bohr:2021din ;
Bjerrum-Bohr:2021vuf ; Bjerrum-Bohr:2021wwt ; Britto:2021pud ;
DiVecchia:2021bdo ; DiVecchia:2021ndb ; Dlapa:2021npj ; Dlapa:2021vgp ;
Herrmann:2021lqe ; Herrmann:2021tct ; Bellazzini:2022wzv ; Bern:2022jvn ;
Blumlein:2022qjy ; DiVecchia:2022piu ; Dlapa:2022lmu ; Heissenberg:2022tsn ;
Jones:2022aji ; Barack:2023oqp ; Bern:2023ccb ; Bjerrum-Bohr:2023iey ;
Damgaard:2023ttc ; Jones:2023ugm ; Kosmopoulos:2023bwc ; Rettegno:2023ghr ;
Wilson-Gerow:2023syq ; Porto:2005ac ; Vaidya:2014kza ; Bini:2017xzy ;
Cachazo:2017jef ; Vines:2017hyw ; Guevara:2017csg ; Bini:2018ywr ;
Guevara:2018wpp ; Levi:2018nxp ; Vines:2018gqi ; Chung:2019duq ;
Damgaard:2019lfh ; Guevara:2019fsj ; Maybee:2019jus ; Aoude:2020onz ;
Bern:2020buy ; Guevara:2020xjx ; Levi:2020kvb ; Levi:2020uwu ;
Jakobsen:2022fcj ; Aoude:2023dui ; Aoude:2023vdk ; DiVecchia:2023frv ;
Heissenberg:2023uvo ; Jakobsen:2023ndj ; Jakobsen:2023hig  along with the
worldline QFT method Bel:1981be ; Westpfahl:1985tsl ; Goldberger:2004jt ;
Goldberger:2006bd ; Ledvinka:2008tk ; Goldberger:2009qd ; Steinhoff:2014kwa ;
Damour:2016gwp ; Blanchet:2018yvb ; Comberiati:2022cpm ; Jakobsen:2022psy ;
Mogull:2020sak ; Saketh:2022xjb ; Wang:2022ntx ; Bhattacharyya:2023kbh ;
Bern:2023ity ; Damgaard:2023vnx ; Kopp:2023eas  has emerged as an efficient
and powerful computational tool that provides a systematic procedure for
converting loop corrections to classical perturbations, which are crucial for
enhancing the accuracy of the description.

However, applying the QFT techniques to binary black hole systems is not
straightforward and poses several challenges. Traditional methods based on the
Feynman diagrams involve cumbersome calculations and do not provide a clear
pathway for incorporating higher-order corrections. Moreover, the perturbative
nature of QFT requires a careful analysis of loop diagrams, which can quickly
become intractable for higher-loop calculations. Modern methods, primarily
based on the double copy Kawai:1985xq ; Bern:2008qj ; BjerrumBohr:2009rd ;
Stieberger:2009hq ; Bjerrum-Bohr:2010diw ; Bjerrum-Bohr:2010mtb ; Bern:2010ue
; Bern:2010yg ; Tye:2010dd ; BjerrumBohr:2010zs ; Feng:2010my ; Mafra:2011kj ;
Monteiro:2011pc ; BjerrumBohr:2012mg ; Monteiro:2014cda ; Luna:2016hge ;
Bern:2019prr ; Bautista:2019tdr  and unitarity cut approach Bern:1994zx ;
Bern:1994cg ; Bern:1995db ; Britto:2004nc ; Ossola:2006us ; Forde:2007mi ;
Berger:2008sj ; Bern:2011qt , offer a more efficient computational framework
that dispenses with explicit Feynman diagrams or Lagrangian.

In the present article, we establish an alternative computational framework
based on the quantum off-shell recursion relations Lee:2022aiu , which is the
loop-level generalization of the Berends-Giele recursions Berends:1987me  for
tree-level. The quantum off-shell recursion provides an efficient and robust
tool for deriving loop integrands. Furthermore, this framework is universally
applicable to any theory with Lagrangian.

We initiate by introducing the field theory of non-rotating binary black hole
systems Goldberger:2004jt ; Porto:2016pyg ; Cheung:2018wkq . In the field
theory description, the non-rotating black holes are realized by massive
scalar fields minimally coupled to gravity. We then construct the perturbative
Lagrangian and its corresponding equations of motion (EoM) by adopting Cheung
and Remmen’s convention Cheung:2016say ; Cheung:2017kzx  of metric
perturbations. This choice yields simpler expressions for the perturbative
expansions than conventional perturbative General Relativity (GR) approaches.
Next, we derive the Dyson-Schwinger (DS) equation, a quantum counterpart to
the classical EoM, by deforming the perturbative EoM. This plays a crucial
role in constructing quantum off-shell recursions.

Subsequently, we construct the quantum perturbiner method Lee:2022aiu  derived
from the quantum effective action formalism. The perturbiner expansion is the
generating function of the off-shell currents directly related to the
scattering amplitude. In this context, the quantum off-shell recursions can be
derived by substituting the perturbiner expansion into the DS equation. We
need to perform the $\hbar$-expansion of off-shell currents, which is crucial
for constructing self-consistent recursion relations. We solve the recursions
up to two-loop order explicitly. Interestingly, the solutions of the off-shell
recursions are loop integrands and include both quantum and classical
contributions.

Next, we introduce a new prescription for extracting and classifying classical
loop diagrams compatible with the off-shell recursion. Since the off-shell
recursion generates loop integrands by solving the recursion relations without
considering explicit diagrams, we need to invent a new prescription, a more
algebraic method. We developed a power-counting prescription by defining an
identifier for diagrams that can be easily implemented in the loop integrands.
To this end, we introduce a set of coupling constants by dividing the Newton
constant $\kappa$ according to the interaction types and a label for each
propagator. We then define the identifier as power counting of the new
couplings and labels of a diagram. It is straightforward to filter out the
integrands corresponding to the classical diagrams from the total result.

Further, we devise a classification scheme for optimizing the IBP reduction
Chetyrkin:1981qh ; Kotikov:1990kg ; Kotikov:1991pm ; Remiddi:1997ny ;
Henn:2013pwa ; Henn:2014qga  process. We first group the classical diagrams
according to the power countings of the new coupling constants, and we call
the groups _classes_. We next categorize each class according to the diagram
topologies. We first define the _primary diagrams_ containing the maximum
number of internal lines. We then may generate diagrams by pinching the
internal lines from a primary diagram. The set of diagrams generated from the
same primary diagram is called _sector_. The main advantage of this
classification is that the diagrams included in the same class incorporate the
same integral family. In the one-loop case, there is only one class with two
sectors. In the two-loop case, there are four classes, and these classes have
their own sectors.

Even for the classical diagrams, numerators of their loop integrands may
include terms that are irrelevant to the classical long-range interaction.
Utilizing the soft expansion that scales the loop momenta and momentum
transfer Parra-Martinez:2020dzs , we truncate the unnecessary terms in the
loop integrands by the power counting rules, thereby obtaining the classical
loop integrands. We also define the soft integral family by replacing the
conventional scalar propagators with linear ones and derive the master
integrals for each sector using IBP reduction. Using the IBP reduction once
more, we may express the classical loop integrands for each sector in terms of
these master integrals, which can be computed by solving the associated
differential equations. We refer to the solutions of the differential
equations derived in Parra-Martinez:2020dzs . Finally, by summing
contributions from each sector, we acquire the explicit expression of the
classical scalar $2\to 2$ amplitudes up to 3PM order. These results are the
same as the known results.

This article is structured as follows: Section 2 introduces our field theory
and derives the perturbative EoMs. We also construct the DS equations using
the deformation of the fields to operators. Section 3 introduces the
perturbiner expansion at the loop level and the quantum off-shell recursion
relations. Section 4 presents our power-counting prescription for filtering
and classifying classical loop diagrams. Section 5 discusses the soft
expansion and IBP reduction techniques with the one-loop example. Section 6
provides explicit computations of scattering amplitudes at two loop levels,
and Section 7 concludes the paper with a summary and outlook for future
research.

##  2 Field Theory for Non-spinning Binary Black Holes

In this section, we discuss the field theory description of a non-rotating
black hole merger in the inspiral phase. During this phase, black holes are
sufficiently separated relative to the horizon scale, thereby allowing us to
approximate the black holes as weakly interacting massive point particles. We
first introduce the field theory, which consists of the Einstein-Hilbert (EH)
action and the minimally coupled massive scalar fields. Next, we derive the
perturbations of the field theory action and its equations of motion (EoM).
Finally, we construct the Dyson-Schwinger (DS) equation, which is crucial in
constructing the quantum off-shell recursion.

###  2.1 Perturbations of the action and EoM

We consider the QFT framework of the non-rotating binary black hole merger in
the inspiral phase. The field theory action consists of two parts: the EH
action $S_{\rm EH}$ and the massive scalar fields $S_{\rm matter}$

$S_{\rm total}=S_{\rm EH}+S_{\rm matter}\,.$ (1)

The dynamics of gravitons are governed by the EH action, which is proportional
to the Ricci scalar $R$

$\displaystyle S_{\rm EH}$ $\displaystyle=\int\mathrm{d}^{4}xL_{\rm
EH}\,,\qquad L_{\rm EH}=\frac{1}{2\kappa^{2}}\sqrt{-g}R\,,$ (2)

where $\kappa^{2}=8\pi G$ and $G$ is the Newton constant. The matter action is
given by massive scalar fields, $\phi_{1}$ and $\phi_{2}$, interacting through
only the minimally coupled gravitons. We also ignore the self-interaction
terms due to the absence of pair-creations or pair-annihilations of black
holes,

$S_{\rm
matter}=\int_{x}\sqrt{-g}\sum_{I=1}^{2}\left(-\frac{1}{2}g^{\mu\nu}\partial_{\mu}\phi_{I}\partial_{\nu}\phi_{I}-\frac{1}{2}m_{I}^{2}\phi^{2}_{I}\right)\,.$
(3)

For a concise perturbation formulation of GR, we transform the metric tensor
to the tensor density in the Landau-Lifshitz book Landau:1975pou . The tensor
density $\sigma$ is defined as

$\sigma^{\mu\nu}=\sqrt{-g}g^{\mu\nu}\,,\qquad\sigma_{\mu\nu}=\frac{1}{\sqrt{-g}}g_{\mu\nu}\,,$
(4)

and $\sigma^{\mu\nu}$ is the inverse of $\sigma_{\mu\nu}$. We recast the EH
action and the matter action in terms of the tensor density up to the boundary
term

$\displaystyle S_{\rm EH}$
$\displaystyle=\frac{1}{2\kappa^{2}}\int\mathrm{d}^{D}x\left[\
\frac{1}{4}\sigma^{\mu\nu}\partial_{\mu}\sigma^{\rho\sigma}\partial_{\nu}\sigma_{\rho\sigma}-\frac{1}{2}\sigma^{\mu\nu}\partial_{\mu}\sigma^{\rho\sigma}\partial_{\rho}\sigma_{\nu\sigma}+(D-2)\sigma^{\mu\nu}\partial_{\mu}\hat{d}\partial_{\nu}\hat{d}\
\right]\,,$ (5) $\displaystyle S_{\rm matter}$
$\displaystyle=\int\mathrm{d}^{D}x\sum_{I=1}^{2}\left[-\frac{1}{2}\sigma^{\mu\nu}\partial_{\mu}\phi_{I}\partial_{\nu}\phi_{I}-\frac{1}{2}m_{I}^{2}|\operatorname{det}\sigma|^{-\frac{1}{2}}\phi_{I}^{2}\right]\,,$

where
$\partial_{\mu}\hat{d}=-\frac{1}{4}\sigma^{\rho\sigma}\partial_{\mu}\sigma_{\rho\sigma}$.

Let us consider the perturbations of the field theory around a flat
background. Although we could consider perturbations of the metric tensor in
the usual way, we employ Cheung and Remmen’s convention of metric
perturbations Cheung:2016say ; Cheung:2017kzx ; Cho:2022faq . This convention
yields much simpler expressions than the usual metric perturbation, defined as
a linear fluctuation of $g_{\mu\nu}$ around a trivial background
$\eta_{\mu\nu}$. In contrast, we define metric perturbations as a linear
fluctuation of the inverse tensor density $\sigma^{\mu\nu}$,

$\sigma^{\mu\nu}=\eta^{\mu\nu}-\kappa
h^{\mu\nu}\,,\qquad\sigma_{\mu\nu}=\eta_{\mu\nu}+\sum_{n=1}^{\infty}\kappa^{n}(h^{n})_{\mu\nu}\,,\qquad|h^{\mu\nu}|\ll
1\,.$ (6)

Here the indices are raised or lowered by the background metric $\eta$

Upon this convention, the all-order expression of the perturbative EH action
(5) and the Einstein equation in $\kappa$ are obtained in Cho:2022faq , and
the perturbations of the Lagrangian action is given by

$\displaystyle L_{\rm EH}=$
$\displaystyle\frac{1}{2\kappa^{2}}\sum_{n=0}^{\infty}\bigg{[}-\frac{1}{4}\kappa^{n+2}(\eta-\kappa
h)^{\mu\nu}\sum_{m=0}^{n}{\operatorname{tr}}\big{[}h^{m}\partial_{\mu}hh^{n-m}\partial_{\nu}h\big{]}+\frac{1}{2}\kappa^{n+2}\big{[}\partial_{\nu}hh^{n}\partial_{\mu}h\big{]}^{\mu\nu}$
(7) $\displaystyle\qquad\qquad+\frac{1}{4(D-2)}\kappa^{n+2}(\eta-\kappa
h)^{\mu\nu}\sum_{m=0}^{n}{\operatorname{tr}}\big{[}h^{m}\partial_{\mu}h\big{]}{\operatorname{tr}}\big{[}h^{n-m}\partial_{\nu}h\big{]}\bigg{]}\,.$

We impose the de Donder gauge choice, $\partial_{\mu}h^{\mu\nu}=0$, and the
gauge fixing term is

$S_{\rm
gf}=\int\mathrm{d}^{4}{x}\Big{(}-\frac{1}{4}\partial_{\mu}h^{\mu\rho}\partial^{\nu}h_{\nu\rho}\Big{)}\,.$
(8)

For general two-loop scattering amplitude, including both quantum and
classical contributions, we should expand the EH action to $\kappa^{3}$.
However, we are considering the “classical” two-loop amplitude – it will be
clear the meaning of classical in Section 4 – and it is sufficient to expand
the action up to $\kappa^{2}$-order terms or four-point vertices,

$L_{\rm EH}+L_{\rm
gf}=\int\mathrm{d}^{4}x\left[\frac{1}{2}L^{h}_{2}+\frac{1}{2}\kappa
L^{h}_{3}+\frac{1}{2}\kappa^{2}L^{h}_{4}+j^{\mu\nu}h_{\mu\nu}\right]\,,$ (9)

where we added the bare external source for graviton $j^{\mu\nu}$. Here
$L^{h}_{2}$ is the graviton kinetic term in the presence of the gauge fixing
term (8)

$L^{h}_{2}=-\frac{1}{4}\partial^{\rho}h^{\mu\nu}\partial_{\rho}h_{\mu\nu}+\frac{1}{8}\partial^{\mu}h\partial_{\mu}h\,,$
(10)

where $h=h^{\mu}{}_{\mu}$. The graviton self-interaction terms $L^{h}_{3}$ and
$L^{h}_{4}$ are given by

$\displaystyle L^{h}_{3}$
$\displaystyle=\frac{1}{2}\big{[}\partial_{\nu}hh\partial_{\mu}h\big{]}^{\mu\nu}-\frac{1}{2}{\operatorname{tr}}\big{[}\partial^{\mu}hh\partial_{\mu}h\big{]}+\frac{1}{4}{\operatorname{tr}}\big{[}h\partial_{\mu}h\big{]}\partial^{\mu}h+\frac{1}{4}h^{\mu\nu}{\operatorname{tr}}\big{[}\partial_{\mu}h\partial_{\nu}h\big{]}$
(11)
$\displaystyle\quad-\frac{1}{8}h^{\mu\nu}\partial_{\mu}h\partial_{\nu}h\,,$
$\displaystyle L^{h}_{4}$
$\displaystyle=-\frac{1}{2}{\operatorname{tr}}\big{[}h^{2}\partial^{\mu}h\partial_{\mu}h\big{]}-\frac{1}{4}{\operatorname{tr}}\big{[}h\partial^{\mu}hh\partial_{\mu}h\big{]}+\frac{1}{4}{\operatorname{tr}}\big{[}h^{2}\partial^{\mu}h\big{]}\partial_{\mu}h+\frac{1}{8}{\operatorname{tr}}\big{[}h\partial^{\mu}h\big{]}{\operatorname{tr}}\big{[}h\partial_{\mu}h\big{]}$
$\displaystyle\quad+\frac{1}{2}h^{\mu\nu}{\operatorname{tr}}\big{[}h\partial_{\mu}h\partial_{\nu}h\big{]}-\frac{1}{4}h^{\mu\nu}{\operatorname{tr}}\big{[}h\partial_{\mu}h\big{]}\partial_{\nu}h+\frac{1}{2}\big{[}\partial_{\nu}hh^{2}\partial_{\mu}h\big{]}^{\mu\nu}\,.$

Here we have employed the matrix notation to avoid cumbersome dummy
indices.111 For arbitrary matrices $A^{\mu}{}_{\nu}$, $B^{\mu}{}_{\nu}$ and
$C^{\mu}{}_{\nu}$, we denote their matrix product as
$A^{\mu}{}_{\rho}B^{\rho}{}_{\nu}=\big{[}AB\big{]}^{\mu}{}_{\nu}\qquad
A^{\mu}{}_{\rho}B^{\rho}{}_{\sigma}C^{\sigma}{}_{\nu}=\big{[}ABC\big{]}^{\mu}{}_{\nu}\,,$
(12) and the trace of matrices as
$A^{\mu}{}_{\rho}B^{\rho}{}_{\mu}:={\operatorname{tr}}\big{[}AB\big{]}\,,\qquad
A^{\mu}{}_{\rho}B^{\rho}{}_{\sigma}C^{\sigma}{}_{\mu}:={\operatorname{tr}}\big{[}ABC\big{]}\,.$
(13) One remark is that we do not need to introduce the ghost fields because
the classical diagrams do not allow the graviton loops.

We may rewrite the kinetic term $L^{h}_{2}$ using the kinetic operator
$K^{\mu\nu,\rho\sigma}[x,y]$

$S^{h}_{2}=\int\mathrm{d}^{4}x\frac{1}{2}L^{h}_{2}=\int\mathrm{d}^{4}x\mathrm{d}^{4}y\bigg{[}-\frac{1}{2}h_{\mu\nu}(x)K^{\mu\nu,\rho\sigma}[x,y]h_{\rho\sigma}(y)\bigg{]}\,,$
(14)

where

$K^{\mu\nu,\rho\sigma}[x,y]=-\frac{1}{4}\bigg{(}\eta^{(\mu|\rho|}\eta^{\nu)\sigma}-\frac{1}{2}\eta^{\mu\nu}\eta^{\rho\sigma}\bigg{)}\Box_{y}\delta^{D}(x-y)\,.$
(15)

Then the graviton propagator $D^{\mu\nu,\rho\sigma}[x,y]$ is the inverse of
$K^{\mu\nu,\rho\sigma}[x,y]$

$D^{\mu\nu,\rho\sigma}[x,y]=2P^{\mu\nu,\rho\sigma}D[x,y]\,,$ (16)

where

$D[x,y]=\int\frac{\mathrm{d}^{4}p}{(2\pi)^{4}}\frac{1}{p^{2}}e^{ip\cdot
x}\,,\qquad
P^{\mu\nu,\rho\sigma}=\eta^{(\mu|\rho|}\eta^{\nu)\sigma}-\frac{1}{2}\eta^{\mu\nu}\eta^{\rho\sigma}\,.$
(17)

Let us consider the perturbative action of the minimally coupled massive
scalars. We should keep the terms up to the $\kappa^{3}$-orders, which is one
order higher than the EH action,

$\displaystyle L_{\rm matter}$
$\displaystyle=\sum_{I=1}^{2}\left[-\frac{1}{2}\sigma^{\mu\nu}\partial_{\mu}\phi_{I}\partial_{\nu}\phi_{I}-\frac{1}{2}m_{I}^{2}|\operatorname{det}\sigma|^{-\frac{1}{2}}\phi_{I}^{2}\right]$
(18)
$\displaystyle=\sum_{I=1}^{2}\Bigg{[}-\frac{1}{2}\big{(}\eta^{\mu\nu}-\kappa
h^{\mu\nu}\big{)}\partial_{\mu}\phi_{I}\partial_{\nu}\phi_{I}-\frac{m_{I}^{2}}{2}\bigg{(}1-\frac{\kappa
h}{2}+\frac{\kappa^{2}h^{2}}{8}-\frac{\kappa^{2}h_{\mu\nu}h^{\mu\nu}}{4}\bigg{)}\phi_{I}^{2}$
$\displaystyle\qquad\quad+\kappa^{3}(\cdots)\Bigg{]}\,.$

For later convenience, we introduce the capital Roman indices $M,N,P\cdots$ to
denote a pair of symmetric indices, such as

$h^{M}=h^{\mu\nu}\,,\qquad P^{M,N}=P^{\mu\nu,\rho\sigma}\,,$ (19)

and in general

$v^{MNP\cdots}=v^{(\mu\nu)(\rho\sigma)(\kappa\lambda)\cdots}\,.$ (20)

To avoid the complicated structures of contractions of dummy indices, we
define tensors that consist of $\eta^{\mu\nu}$ and $\delta_{\mu}{}^{\nu}$. We
will call them as $\mathcal{I}$-tensors. For instance, let us consider a
matrix product of three $h^{\mu\nu}$,

$\displaystyle{}[hhh]^{\mu\nu}$
$\displaystyle=h^{\mu\rho}h^{\rho\sigma}h^{\sigma\nu}$ (21)
$\displaystyle=\big{(}\eta^{\mu\kappa_{1}}\eta^{\nu\lambda_{3}}\eta^{\lambda_{1}\kappa_{2}}\eta^{\lambda_{2}\kappa_{3}}\big{)}h^{\kappa_{1}\lambda_{1}}h^{\kappa_{2}\lambda_{2}}h^{\kappa_{3}\lambda_{3}}\,,$
$\displaystyle=\mathcal{I}^{\mu\nu,\kappa_{1}\lambda_{1},\kappa_{2}\lambda_{2},\kappa_{3}\lambda_{3}}h^{\kappa_{1}\lambda_{1}}h^{\kappa_{2}\lambda_{2}}h^{\kappa_{3}\lambda_{3}}\,,$
$\displaystyle=\mathcal{I}^{\mu\nu,N,P,Q}h^{N}h^{P}h^{Q}\,.$

where
$\mathcal{I}^{\mu\nu,\kappa_{1}\lambda_{1},\kappa_{2}\lambda_{2},\kappa_{3}\lambda_{3}}=\eta^{\mu(\kappa_{1}}\eta^{\lambda_{1})(\kappa_{2}}\eta^{\lambda_{2})(\kappa_{3}}\eta^{\lambda_{3})\nu}$.
Hereafter, we will ignore the position of the indices while keeping the
summation convention for repeated indices. In Appendix A, we present the
explicit expressions of the $\mathcal{I}$-tensors and the general structure of
the EoMs.

Equipped with the $\mathcal{I}$-tensors, we now represent the EoMs of
$h^{\mu\nu}$ and $\phi^{I}$. The EoM of $h^{\mu\nu}$ is divided into two
parts: the perturbations of the curvature part $\mathcal{G}^{M}_{n}$ and the
energy-momentum (EM) tensor $\mathcal{T}^{M}_{n}$

$\Box
h^{M}=2P^{M,N}\Biggl{[}\sum_{n=1}^{2}\kappa^{n}\mathcal{G}^{N}_{n}-\sum_{n=1}^{3}\kappa^{n}\mathcal{T}^{N}_{n}-2j^{N}\Bigg{]}\,.$
(22)

Here, curvature perturbations $\mathcal{G}^{M}_{1}$ and $\mathcal{G}^{M}_{2}$
are obtained by variations of $L^{h}_{3}$ and $L^{h}_{4}$ in (11) respectively

$\displaystyle\mathcal{G}^{M}_{1}$ $\displaystyle=\mathcal{I}^{M,\rho N,\sigma
P}_{1,1}\partial_{\rho}h^{N}\partial_{\sigma}h^{P}+\mathcal{I}^{M,N,\rho\sigma
P}_{1,2}h^{N}\partial_{\rho}\partial_{\sigma}h^{P}\,,$ (23)
$\displaystyle\mathcal{G}^{M}_{2}$ $\displaystyle=\mathcal{I}^{M,N,\rho
P,\sigma
Q}_{2,1}h^{N}\partial_{\rho}h^{P}\partial_{\sigma}h^{Q}+\mathcal{I}^{M,N,P,\rho\sigma
Q}_{2,2}h^{N}h^{P}\partial_{\rho}\partial_{\sigma}h^{Q}\,,$

where $\mathcal{I}_{m,n}$ are defined in (232) and (236). The EM tensor parts,
$\mathcal{T}^{\mu\nu}_{1}$, $\mathcal{T}^{\mu\nu}_{2}$ and
$\mathcal{T}^{\mu\nu}_{3}$, are given by variations of the perturbed matter
action (18),

$\displaystyle\mathcal{T}^{\mu\nu}_{1}$
$\displaystyle=\sum_{I=1}^{2}\left[\partial^{\mu}\phi_{I}\partial^{\nu}\phi_{I}+\frac{m_{I}^{2}}{2}\eta^{\mu\nu}\phi_{I}\phi_{I}\right]\,,$
(24) $\displaystyle\mathcal{T}^{M}_{2}$ $\displaystyle=\mathcal{I}_{\rm
EM}^{M,N}\sum_{I=1}^{2}m^{2}_{I}h^{N}\phi_{I}\phi_{I}\,,\qquad\mathcal{T}^{M}_{3}=\mathcal{I}_{\rm
EM}^{M,N,P}\sum_{I=1}^{2}m^{2}_{I}h^{N}h^{P}\phi_{I}\phi_{I}\,,$

where we present $\mathcal{I}_{{\rm EM}}$ in (240) and (241).

Similarly, perturbations of EoM of the scalar fields can be represented as

$\displaystyle\left(\square-m_{I}^{2}\right)\phi_{I}$
$\displaystyle=-\rho_{I}+\kappa\bigg{[}\partial_{\kappa}\big{(}h^{\kappa\lambda}\partial_{\lambda}\phi_{I}\big{)}-\frac{\kappa
m_{I}^{2}}{2}h^{\kappa\kappa}\phi_{I}\bigg{]}$ (25)
$\displaystyle\quad+\kappa^{2}m^{2}_{I}\mathcal{I}^{M,N}_{\phi}h^{M}h^{N}\phi_{I}+\kappa^{3}m^{2}_{I}\mathcal{I}^{M,N,P}_{\phi}h^{M}h^{N}h^{P}\phi_{I}\,,$

where $\mathcal{I}_{\phi}$ are defined in (242).

###  2.2 Dyson-Schwinger equations

The Dyson-Schwinger (DS) equation, a quantum counterpart of the classical EoM,
is crucial in constructing the quantum off-shell recursions. It can be
obtained simply by deforming the fields of the classical EoMs. In this
context, the classical fields are defined by the vacuum expectation value of
the fields

$\displaystyle\varphi^{I}_{x}=\frac{\delta
W[j,\rho]}{\delta\rho^{I}_{x}}\,,\qquad h^{\mu\nu}_{x}=\frac{\delta
W[j,\rho]}{\delta j^{\mu\nu}_{x}}\,,$ (26)

where $\rho^{I}_{x}$ and $j^{\mu\nu}_{x}$ are the bare external sources for
$\phi^{I}$ and $h^{\mu\nu}$, and $W[j,\rho]$ is the generating functional for
connected diagrams

$\exp\big{[}W[j,\rho]\big{]}=\int\mathcal{D}\phi^{I}\mathcal{D}h_{\mu\nu}\exp\bigg{[}\frac{i}{\hbar}\Big{(}S_{\rm
EH}+S_{\rm
matter}+\rho^{I}_{x}\phi^{I}_{x}+j^{\mu\nu}_{x}h_{x}^{\mu\nu}\Big{)}\bigg{]}\,.$
(27)

We will encounter the dressed external sources in the next section. Here, the
subscripts denote coordinates such as

$\varphi^{I}_{x}=\varphi^{I}(x)\,,\qquad\rho^{I}_{x}=\rho^{I}(x)\,,\qquad
h^{\mu\nu}_{x}=h^{\mu\nu}(x)\,,\qquad j^{\mu\nu}_{x}=j^{\mu\nu}(x)\,.$ (28)

We will adopt this notation throughout this paper.

The DS equations are generated by the following deformations of the fields
from the classical EoM Lee:2022aiu

$\displaystyle\phi^{I}_{x}\to\varphi^{I}+\frac{\hbar}{i}\frac{\delta}{\delta\rho^{I}_{x}}\,,\qquad
h^{M}_{x}\to h^{\mu\nu}_{x}+\frac{\hbar}{i}\frac{\delta}{\delta
j^{\mu\nu}_{x}}\,.$ (29)

We summarized the general structure of the deformation rule in Appendix B.1.

Thus, the DS equation includes several functional derivatives of the
$\varphi^{I}$ and $h^{\mu\nu}$, which makes it difficult to solve. We follow
the strategy that treats the functional derivatives as independent degrees of
freedom, called _descendant fields_ Lee:2022aiu . Considering up to two-loop
order, the first and the second descendant fields arise in the DS equations.
The first descendant fields are the single functional derivatives of the
classical fields

$\displaystyle\psi^{MN}_{xy}$ $\displaystyle=\frac{\delta h_{x}^{M}}{\delta
j^{N}_{y}}\,,$ $\displaystyle\quad\psi^{MI}_{xy}$ $\displaystyle=\frac{\delta
h_{x}^{M}}{\delta\rho^{I}_{y}}\,,$ $\displaystyle\quad\psi^{IM}_{xy}$
$\displaystyle=\frac{\delta\varphi_{x}^{I}}{\delta j^{M}_{y}}\,,$
$\displaystyle\quad\psi^{IJ}_{xy}$
$\displaystyle=\frac{\delta\varphi_{x}^{I}}{\delta\rho^{J}_{y}}\,,$ (30)

and the second descendant fields are the double-functional derivatives of the
classical fields

$\displaystyle\psi^{MNP}_{xyz}$
$\displaystyle=\frac{\delta^{2}h_{x}^{M}}{\delta j^{N}_{y}\delta
j^{P}_{z}}\,,$ $\displaystyle\quad\psi^{MNI}_{xyz}$
$\displaystyle=\frac{\delta^{2}h_{x}^{M}}{\delta
j^{N}_{y}\delta\rho^{I}_{z}}\,,$ $\displaystyle\quad\psi^{MIN}_{xyz}$
$\displaystyle=\frac{\delta^{2}h_{x}^{M}}{\delta\rho^{I}_{y}\delta
j^{N}_{z}}\,,$ $\displaystyle\quad\psi^{MIJ}_{xyz}$
$\displaystyle=\frac{\delta^{2}h_{x}^{M}}{\delta\rho^{I}_{y}\delta\rho^{J}_{z}}\,,$
(31) $\displaystyle\psi^{IMN}_{xyz}$
$\displaystyle=\frac{\delta^{2}\varphi_{x}^{I}}{\delta j^{M}_{y}\delta
j^{N}_{z}}\,,$ $\displaystyle\quad\psi^{IMJ}_{xyz}$
$\displaystyle=\frac{\delta^{2}\varphi_{x}^{I}}{\delta
j^{M}_{y}\delta\rho^{J}_{z}}\,,$ $\displaystyle\quad\psi^{IJM}_{xyz}$
$\displaystyle=\frac{\delta^{2}\varphi_{x}^{I}}{\delta\rho^{J}_{y}\delta
j^{M}_{z}}\,,$ $\displaystyle\quad\psi^{IJK}_{xyz}$
$\displaystyle=\frac{\delta^{2}\varphi_{x}^{I}}{\delta\rho^{J}_{y}\rho^{K}_{z}}\,.$

Again, for concise expressions, we introduce collective indices $A,B,C,\cdots$
by combining $M,N,\cdots$ and $I,J,\cdots$,

$A=\\{M,I\\}\,,\quad B=\\{N,J\\}\,,\quad\cdots\,.$ (32)

Using the notations, we may define the collective external sources
$\mathbf{j}^{A}=\\{j^{\mu\nu},\rho^{I}\\}=\\{j^{M},\rho^{I}\\}$. The
collective classical fields are given by

$\mathbf{h}^{A}_{x}=\frac{\delta
W[\mathbf{j}]}{\delta\mathbf{j}^{A}_{x}}=\big{\\{}h^{M}_{x},\varphi^{I}_{x}\big{\\}}\,,$
(33)

and their descendants defined in (30) and (31) are represented in a simple
form

$\displaystyle\psi^{IA}_{xy}=\frac{\delta\varphi_{x}^{I}}{\delta\mathbf{j}_{y}^{A}}\,,\qquad\psi^{MA}_{xy}=\frac{\delta
h^{M}_{x}}{\delta\mathbf{j}_{y}^{A}}\,,\qquad\psi^{IAB}_{xyz}=\frac{\delta^{2}\varphi_{x}^{I}}{\delta\mathbf{j}^{A}_{y}\delta\mathbf{j}^{B}_{z}}\,,\qquad\psi^{MAB}_{xyz}=\frac{\delta
h^{M}_{x}}{\delta\mathbf{j}^{A}_{y}\delta\mathbf{j}^{B}_{z}}\,.$ (34)

In solving the DS equations, the $\hbar$-expansion is crucial. Since the
descendant fields are treated by independent degrees of freedom, their field
equations should be imposed for solving the DS equation. We may derive the EoM
for the descendant fields by acting on additional functional derivatives to
the DS equation. However, this process never stops, and one encounters an
infinite sequence of new descendant fields upon each functional derivative. To
bypass this issue, we should expand the DS equation to a particular loop
order. Specifically, we can expand the DS equation in $\hbar$, directly
related to the loop order. Notably, all the classical fields and their
descendant fields are derived from the full generating functional
$W[\mathbf{j}]$ without any approximation; these contain all order expansions
in $\hbar$

$\displaystyle\mathbf{h}^{A}_{x}=\sum_{n=0}^{\infty}\mathbf{h}^{{\scriptscriptstyle(n)}}|^{A}_{x}\hbar^{n}\,,\qquad\psi^{A,B}_{x,y}=\sum_{n=0}^{\infty}\psi^{{\scriptscriptstyle(n)}}|^{A,B}_{x,y}\hbar^{n}\,,\qquad\psi^{A,B,C}_{x,y,z}=\sum_{n=0}^{\infty}\psi^{{\scriptscriptstyle(n)}}|^{A,B,C}_{x,y,z}\hbar^{n}\,.$
(35)

If we substitute these $\hbar$-expansions and keep the terms at a fixed order
in $\hbar$, we can truncate the new descendant fields because these are higher
$\hbar$-order terms.

The DS equation for the scalar fields is represented by the
$\mathcal{I}_{\phi}$ tensors as

$\displaystyle\left(\square-
m_{I}^{2}\right)\varphi^{I}_{x}=-\rho_{x}^{I}+\kappa\Omega^{I}_{1}+\kappa^{2}\mathcal{I}^{I,M,N}_{\phi}\Omega^{IMN}_{2}+\kappa^{3}\mathcal{I}^{I,M,N,P}_{\phi}\Omega^{IMNP}_{3}\,,$
(36)

where the explicit form of the $\Omega$-tensors are listed in (248). It is
straightforward to derive the field equations for the descendant fields by
acting functional derivative $\frac{\delta}{\delta\mathbf{j}^{A}_{x}}$ on the
DS equation (36). One can find their explicit form in (255) and (258).

The DS equation for $h^{M}$ is given by

$\Box
h_{x}^{M}=2P^{M,N}\bigg{(}\sum_{n=1}^{2}\kappa^{n}\mathcal{G}^{M}_{n}-\sum_{n=1}^{3}\kappa^{n}\mathcal{T}^{N}_{n}-2j^{N}_{x}\bigg{)}\,,$
(37)

where $\mathcal{G}^{M}_{x}$ represents the curvature perturbations

$\displaystyle\mathcal{G}^{M}_{1}$ $\displaystyle=\mathcal{I}^{M,\rho N,\sigma
P}_{1,1}\Lambda^{\rho N,\sigma P}_{1,1}+\mathcal{I}^{M,N,\rho\sigma
P}_{1,2}\Lambda^{N,\rho\sigma P}_{1,2}\,,$ (38)
$\displaystyle\mathcal{G}^{M}_{2}$ $\displaystyle=\mathcal{I}^{M,N,\rho
P,\sigma Q}_{2,1}\Lambda^{N,\rho P,\sigma
Q}_{2,1}+\mathcal{I}^{M,N,P,\rho\sigma Q}_{2,2}\Lambda^{N,P,\rho\sigma
Q}_{2,2}\,,$

and

$\displaystyle\mathcal{T}^{M}_{1}=\Theta^{M}_{1}\,,\qquad\mathcal{T}^{M}_{2}=\mathcal{I}_{\rm
EM}^{M,N}\Theta^{N}_{2}\,,\qquad\mathcal{T}^{M}_{3}=\mathcal{I}_{\rm
EM}^{M,N,P}\Theta^{N,P}_{3}\,.$ (39)

The explicit form of the $\Lambda$ and $\Theta$ fields are listed in (253) and
(254).

As we have stated, the field equations for the descendant fields are generated
through the successive applications of the functional derivatives
$\frac{\delta}{\delta\mathbf{j}^{A}_{x}}$ to the DS equations. The
comprehensive formulations of these equations are presented in (260) and
(265).

##  3 Quantum Perturbiner methods for Binary Black Hole System

In the present section, we will derive the quantum off-shell recursion
relations up to the two-loop level. The central object in the off-shell
recursion is the off-shell current, which is given by a correlation function
with at least an off-shell leg and amputated on-shell external legs. We first
construct the quantum perturbiner expansion for $\mathbf{h}^{A}_{x}$ and its
descendant fields, employing the quantum effective action formalism. We then
derive the recursions by substituting these expansions into the DS equations
derived in the previous section. We also discuss the initial conditions of the
recursions.

###  3.1 Quantum perturbiner expansion

The perturbiner method was introduced for deriving the off-shell recursions at
tree level using the classical EoM Selivanov:1997an ; Selivanov:1997aq ;
Rosly:1996vr ; Rosly:1997ap ; Selivanov:1997ts  instead of patterns in
interaction vertices. It has also been applied to gravity Gomez:2021shh ;
Cho:2022faq ; Tao:2023yxy . Recently, the perturbiner method has been
generalized to the quantum perturbiner method by including loop corrections
Lee:2022aiu ; Gomez:2022dzk ; Chen:2023bji . Here we will employ the approach
based on the DS equation Lee:2022aiu , which is available for higher loops. In
this subsection, we derive the quantum perturbiner expansion from the
definition of the classical fields (33) and their descendants.

Let us consider the formal series expansion of the generating functional for
connected diagrams $W[\mathbf{j}]$ respective to the collective dressed
external source $\mathbf{j}^{A}_{x}$

$\displaystyle W[\mathbf{j}]$
$\displaystyle=\sum_{n>1}\frac{1}{n!}\int_{x_{1}x_{2}\cdots
x_{n}}\frac{\delta^{n}W[\mathbf{j}]}{\delta\mathbf{j}^{A_{1}}_{x_{1}}\delta\mathbf{j}^{A_{2}}_{x_{2}}\cdots\delta\mathbf{j}^{A_{n}}_{x_{n}}}\Big{|}_{\mathbf{j}=0}\mathbf{j}^{A_{1}}_{x_{1}}\mathbf{j}^{A_{2}}_{x_{2}}\cdots\mathbf{j}^{A_{n}}_{x_{n}}\,,$
(40)
$\displaystyle=\frac{\hbar}{i}\sum_{n>1}\frac{1}{n!}\int_{x_{1}x_{2}\cdots
x_{n}}\langle
0\\!\left|T\big{[}\mathbf{h}^{A_{1}}_{x_{1}}\mathbf{h}^{A_{2}}_{x_{2}}\cdots\mathbf{h}^{A_{n}}_{x_{n}}\big{]}\right|\\!0\rangle_{c}\frac{i\mathbf{j}^{A_{1}}_{x_{1}}}{\hbar}\frac{i\mathbf{j}^{A_{2}}_{x_{2}}}{\hbar}\cdots\frac{i\mathbf{j}^{A_{n}}_{x_{n}}}{\hbar}\,,$

where $\langle
0\\!\left|T\big{[}\mathbf{h}^{A_{1}}_{x_{1}}\mathbf{h}^{A_{2}}_{x_{2}}\cdots\mathbf{h}^{A_{n}}_{x_{n}}\big{]}\right|\\!0\rangle_{c}$
is a connected time-ordered $n$-point correlator and

$\int_{x,y\cdots}=\int\mathrm{d}^{4}x\mathrm{~{}d}^{4}y\cdots\quad\text{ and
}\quad\int_{p,q,\cdots}=\int\frac{\mathrm{d}^{4}p}{(2\pi)^{4}}\frac{\mathrm{d}^{4}q}{(2\pi)^{4}}\cdots\,.$
(41)

The corresponding one-point function $\frac{\delta
W[\mathbf{j}]}{\delta\mathbf{j}^{A}_{x}}$ in the presence of the external
source is given by

$\displaystyle\mathbf{h}^{A}_{x}=\sum_{n=1}^{\infty}\frac{1}{n!}\int_{y_{1},y_{2},\cdots,y_{n}}\langle
0|T\big{[}\mathbf{h}^{A}_{x}\mathbf{h}^{B_{1}}_{y_{1}}\cdots\mathbf{h}^{B_{n}}_{y_{n}}\big{]}|0\rangle_{c}\frac{i\mathbf{j}^{B_{1}}_{y_{1}}}{\hbar}\cdots\frac{i\mathbf{j}^{B_{n}}_{y_{n}}}{\hbar}\,.$
(42)

We next identify the external source $\mathbf{j}^{A}_{x}$ to make contact with
the scattering amplitude. According to the LSZ reduction formula,
$W[\mathbf{j}]$ can be written as the summation of scattering amplitudes by
choosing the external source $\mathbf{j}^{A}_{x}$ as follows:

$\mathbf{j}^{A}_{x}=\sum_{i=1}^{N-1}\int_{y_{i}}\mathbf{K}^{AB_{i}}_{xy_{i}}\boldsymbol{\epsilon}^{B_{i}}_{i}e^{-ik_{i}\cdot
y_{i}}=\sum_{i=1}^{N-1}\tilde{\mathbf{K}}^{AB_{i}}_{k_{i}}\boldsymbol{\epsilon}^{B_{i}}_{i}e^{-ik_{i}\cdot
x}\,,$ (43)

where $N$ is the total number of external particles, and
$\mathbf{K}^{AB}_{xy}$ is the inverse of the dressed propagator
$\mathbf{D}^{AB}_{xy}$ in the position space
($\tilde{\mathbf{K}}^{AB_{i}}_{k_{i}}$ and $\tilde{\mathbf{D}}^{AB}_{k_{i}}$
in momentum space) including all-loop corrections, which amputates the
external lines. These satisfy the inverse relation:

$\displaystyle\int_{y}\mathbf{D}^{AC}_{x,y}\mathbf{K}^{CB}_{y,z}=\delta^{AB}\delta_{x,z}\,,\qquad\text{or}\qquad\tilde{\mathbf{K}}^{AC}(p)\tilde{\mathbf{D}}^{CB}(p)=\delta^{AB}\,.$
(44)

Here $\boldsymbol{\epsilon}^{A}_{i}$ is the collective polarization tensor
that consists of the graviton polarization tensor $\epsilon^{\mu\nu}_{i}$ and
scalar polarization $\epsilon^{I}_{i}$

$\boldsymbol{\epsilon}^{A}_{i}=\big{\\{}\epsilon^{\mu\nu}_{i},\epsilon^{I}_{i}\big{\\}}\,.$
(45)

The scalar polarization decides the type of the scalar fields $\varphi^{1}$ or
$\varphi^{2}$ for each external scalar state. The choice of the polarization
depends on the problems in which we are interested. Since we are considering
$2\to 2$ scattering between two scalar particles, we choose the polarizations
as follows:

$\displaystyle\epsilon^{1}_{1}$ $\displaystyle=1\,,$
$\displaystyle\qquad\epsilon^{1}_{2}$ $\displaystyle=0\,,$
$\displaystyle\qquad\epsilon^{1}_{3}$ $\displaystyle=0\,,$ (46)
$\displaystyle\epsilon^{2}_{1}$ $\displaystyle=0\,,$
$\displaystyle\qquad\epsilon^{2}_{2}$ $\displaystyle=1\,,$
$\displaystyle\qquad\epsilon^{2}_{3}$ $\displaystyle=1\,,$
$\displaystyle\epsilon^{\mu\nu}_{1}$ $\displaystyle=0\,,$
$\displaystyle\qquad\epsilon^{\mu\nu}_{2}$ $\displaystyle=0\,,$
$\displaystyle\qquad\epsilon^{\mu\nu}_{3}$ $\displaystyle=0\,.$

This is compatible with the kinematics depicted in Figure 2.

One remark is that we should distinguish the dressed external sources from the
bare ones in the action and DS equations. The bare sources are written in the
inverse propagators or the kinetic operators,

$j^{M}_{x}=\sum_{i=1}^{N-1}\tilde{K}^{MP}_{k_{i}}\epsilon^{P}_{i}e^{-ik_{i}\cdot
x}\,,\qquad\rho^{I}_{x}=\sum_{i=1}^{N-1}\tilde{K}^{IJ}_{k_{i}}\epsilon^{J}_{i}e^{-ik_{i}\cdot
x}\,,$ (47)

where $\tilde{K}^{\mu\nu,\rho\sigma}_{k_{i}}$ and $\tilde{K}^{IJ}_{k_{i}}$ are
the kinetic operators defined by

$\displaystyle\tilde{K}^{\mu\nu,\rho\sigma}_{k_{i}}$
$\displaystyle=\frac{1}{2}\bigg{(}\eta^{\mu\rho}\eta^{\nu\sigma}-\frac{1}{2}\eta^{\mu\nu}\eta^{\rho\sigma}\bigg{)}k_{i}^{2}\,,$
(48) $\displaystyle\tilde{K}^{IJ}_{k_{i}}$
$\displaystyle=\delta^{IJ}\big{(}k_{i}^{2}+m_{I}^{2}\big{)}\,.$

After substituting the definition of the dressed external source in (43) into
the expansion of $\mathbf{h}^{A}_{x}$ (42), we convert to the momentum space
by Fourier transformation

$\displaystyle\mathbf{h}^{A}_{x}=\sum_{n=1}^{\infty}\sum_{i_{1},\cdots
i_{n}=1}^{n}\int_{\ell}$
$\displaystyle(2\pi)^{4}\delta^{4}\big{(}l+k_{i_{1}\cdots i_{n}}\big{)}\langle
0|\tilde{\mathbf{h}}^{A}(l)\tilde{\mathbf{h}}^{B_{i_{1}}}(k_{i_{1}})\cdots\tilde{\mathbf{h}}^{B_{i_{n}}}(k_{i_{n}})|0\rangle_{c}$
(49)
$\displaystyle\times\Big{(}\frac{i}{\hbar}\Big{)}^{n}\big{(}\tilde{\mathbf{K}}^{B_{i_{1}}C_{i_{1}}}_{k_{i_{1}}}\boldsymbol{\epsilon}^{C_{i_{1}}}\big{)}\cdots\big{(}\tilde{\mathbf{K}}^{B_{i_{n}}C_{i_{n}}}_{k_{i_{n}}}\boldsymbol{\epsilon}^{C_{i_{n}}}\big{)}e^{il\cdot
x}\,,$

where $k_{i_{1}\cdots i_{n}}=k_{i_{1}}+k_{i_{2}}+\cdots+k_{i_{n}}$. We may
recast the expansion and define the _quantum perturbiner expansion_
Lee:2022aiu  for $\mathbf{h}^{A}_{x}$

$\displaystyle\mathbf{h}^{A}_{x}$
$\displaystyle=\sum_{n=1}^{\infty}\sum_{i_{1},\cdots,i_{n}=1}^{n}\mathbf{J}^{A}_{i_{1}\cdots
i_{n}}e^{-ik_{i_{1}i_{2}\cdots i_{n}}\cdot x}\,,$ (50)

where the coefficients $\mathbf{J}^{A}_{i_{1}\cdots i_{n}}$ is the _quantum
off-shell currents_ and its explicit form is

$\displaystyle\mathbf{J}^{A}_{i_{1}\cdots
i_{n}}=\big{\langle}0\big{|}\tilde{\mathbf{h}}^{A}(-k_{i_{1}i_{2}\cdots
i_{n}})\tilde{\mathbf{h}}^{B_{1}}(k_{i_{1}})\tilde{\mathbf{h}}^{B_{2}}(k_{i_{2}})\cdots\tilde{\mathbf{h}}^{B_{n}}(k_{i_{n}})\big{|}0\big{\rangle}_{c}$
(51)
$\displaystyle\times\Big{(}\frac{i}{\hbar}\Big{)}^{n}\big{(}\tilde{\mathbf{K}}^{B_{i_{1}}C_{i_{1}}}_{k_{i_{1}}}\boldsymbol{\epsilon}^{C_{i_{1}}}\big{)}\cdots\big{(}\tilde{\mathbf{K}}^{B_{i_{n}}C_{i_{n}}}_{k_{i_{n}}}\boldsymbol{\epsilon}^{C_{i_{n}}}\big{)}\,.$

This implies the quantum perturbiner expansion is a generating function of the
off-shell currents $\mathbf{J}^{A}_{i_{1}\cdots i_{n}}$. Apparently, the
definition of the quantum perturbiner expansion looks the same as the
classical one, but the quantum case involves all order corrections in $\hbar$.

It is convenient to denote the quantum perturbiner expansion (50) in terms of
the so-called _words_.

$\mathbf{h}^{A}_{x}=\sum_{\mathcal{P}}\mathbf{J}^{A}_{\mathcal{P}}e^{-ik_{\mathcal{P}}\cdot
x}\,,$ (52)

where $\mathcal{P}$ is the word labelling multi-particle states and
$\sum_{\mathcal{P}}$ is a summation over all possible words. In the
perturbiner formalism, single-particle states are labeled by what we refer to
as “letters,” denoted by $i,j,k,l,\cdots$, whereas combinations of these
states, or “words,” are represented as
$\mathcal{P},\mathcal{Q},\mathcal{R}\cdots$, for instance, $\mathcal{P}=ijk$.
Notably, within our field theory description, the order of letters within a
given word is meaningless as single-particle states are represented solely by
their momenta $k_{\mathcal{P}}$ in the plane wave basis. For instance,
$\mathbf{J}^{A}_{123}=\mathbf{J}^{A}_{132}=\cdots$. To eliminate such
redundancy, we impose a canonical ordering on the letters within each word,
such that for an arbitrary word $\mathcal{P}=ijkl\cdots$, we have
$i<j<k<l\cdots$.

It is evident that off-shell currents associated with the empty word, denoted
by $\mathcal{P}=\emptyset$, are trivial, leading to

$\mathbf{J}^{A}_{\emptyset}=0\,.$ (53)

On the other hand, as we will elaborate below, descendant currents do admit
zero modes. The length of a word is called its “rank” denoted as
$|\mathcal{P}|,|\mathcal{Q}|$, $|\mathcal{R}|$. We further require that the
off-shell currents with repeated momenta vanish,

$\mathbf{h}^{A}_{i_{1}\cdots j\cdots j\cdots i_{n}}=0\,.$ (54)

Whenever we want, we can separate the perturbiner expansion for the collective
currents $\mathbf{h}^{A}_{x}$ into its components $h^{\mu\nu}_{x}$ and
$\varphi^{I}_{x}$

$\displaystyle
h^{\mu\nu}_{x}=\sum_{\mathcal{P}}J^{\mu\nu}_{\mathcal{P}}e^{-ik_{\mathcal{P}}\cdot
x}\,,\qquad\varphi^{I}_{x}$
$\displaystyle=\sum_{\mathcal{P}}\Phi^{I}_{\mathcal{P}}e^{-ik_{\mathcal{P}}\cdot
x}\,,$ (55)

where $\Phi^{I}_{\mathcal{P}}$ and $J^{\mu\nu}_{\mathcal{P}}$ are the scalars
and graviton off-shell currents respectively.

According to the LSZ reduction formula, the off-shell current is directly
related to the scattering amplitudes after amputating the off-shell leg and
taking the on-shell limit. For instance, the $n$-loop four-point scattering
amplitude between the scalar particles is given by

$\left(\frac{\hbar}{i}\right)^{n}\mathcal{M}^{(n)}\left(k_{1},k_{2},k_{3},k_{4}\right)=\lim_{k_{123}^{2}\rightarrow-m^{2}_{1}}\sum_{p=0}^{n}\tilde{\mathbf{K}}^{(p)}_{123}\Phi^{(n-p)}\big{|}^{1}_{123}\,,$
(56)

where $k_{i}$ are external momenta.

We now consider the quantum perturbiner expansion for the descendant fields.
The definition of the descendant fields can be rewritten by multiple
functional derivatives on $W[\mathbf{j}]$, which is a $n$-point function in
the presence of the external source,

$\psi^{A_{1}A_{2}\cdots A_{n}}_{x_{1}x_{2}\cdots
x_{n}}=\frac{\delta^{n}W[\mathbf{j}]}{\delta\mathbf{j}^{A_{1}}_{x_{1}}\delta\mathbf{j}^{A_{2}}_{x_{2}}\cdots\delta\mathbf{j}^{A_{n}}_{x_{n}}}\,.$
(57)

Similar to (49), descendant fields can be represented in terms of the
correlation function with inverse two-point functions such as

$\displaystyle\psi^{A_{1}A_{2}\cdots A_{n}}_{x_{1}x_{2}\cdots x_{n}}=$
$\displaystyle\sum_{m=1}^{\infty}\sum_{j_{1},\cdots
j_{m}=1}^{m}\int_{\ell_{1}\ell_{2}\cdots\ell_{n}}(2\pi)^{4}\delta^{4}\big{(}\ell_{1}+\ell_{2}+\cdots+\ell_{n}+k_{j_{1}\cdots
j_{m}}\big{)}$ (58) $\displaystyle\times\langle
0|\tilde{\mathbf{h}}^{A_{1}}(\ell_{1})\cdots\tilde{\mathbf{h}}^{A_{n}}(\ell_{n})\tilde{\mathbf{h}}^{B_{j_{1}}}(k_{j_{1}})\cdots\tilde{\mathbf{h}}^{B_{j_{m}}}(k_{j_{m}})|0\rangle_{c}$
$\displaystyle\times\Big{(}\frac{i}{\hbar}\Big{)}^{m}\big{(}\tilde{\mathbf{K}}^{B_{j_{1}}C_{j_{1}}}_{k_{j_{1}}}\boldsymbol{\epsilon}^{C_{j_{1}}}\big{)}\cdots\big{(}\tilde{\mathbf{K}}^{B_{j_{m}}C_{j_{m}}}_{k_{j_{m}}}\boldsymbol{\epsilon}^{C_{j_{m}}}\big{)}e^{i(\ell_{1}\cdot
x_{1}+\cdots+\ell_{n}\cdot x_{n})}\,.$

We define the $n$-th descendant currents $\Psi^{A_{1}A_{2}\cdots
A_{n+1}}_{\ell_{1}\ell_{2}\cdots\ell_{n}|\hat{\mathcal{P}}}$ in terms of the
connected correlation functions

$\displaystyle\Psi^{A_{1}A_{2}\cdots
A_{n+1}}_{\ell_{1}\ell_{2}\cdots\ell_{n}|\hat{\mathcal{P}}}=\sum_{m=1}^{\infty}\sum_{i_{1},\cdots
i_{m}=1}^{m}\int_{\ell_{1}\ell_{2}\cdots\ell_{n}}\Big{(}\frac{i}{\hbar}\Big{)}^{m}\big{(}\tilde{\mathbf{K}}^{B_{i_{1}}C_{i_{1}}}_{k_{i_{1}}}\boldsymbol{\epsilon}^{C_{i_{1}}}\big{)}\cdots\big{(}\tilde{\mathbf{K}}^{B_{i_{m}}C_{i_{m}}}_{k_{i_{m}}}\boldsymbol{\epsilon}^{C_{i_{m}}}\big{)}$
(59) $\displaystyle\quad\times\langle
0|\tilde{\mathbf{h}}^{A_{1}}(\ell_{1})\cdots\tilde{\mathbf{h}}^{A_{n}}(\ell_{n})\tilde{\mathbf{h}}^{A_{n+1}}(-\ell_{1\cdots
n}-k_{1\cdots
n})\tilde{\mathbf{h}}^{A_{n}}(\ell_{n})\tilde{\mathbf{h}}^{B_{i_{1}}}(k_{i_{1}})\cdots\tilde{\mathbf{h}}^{B_{i_{m}}}(k_{i_{m}})|0\rangle_{c}\,,$

where $\hat{\mathcal{P}}$ denotes a word including the emptyset such as
$\hat{\mathcal{P}}=\mathcal{P}\cup\emptyset$. Upon this definition, we derive
the perturbiner expansion

$\psi^{A_{1}A_{2}\cdots A_{n}}_{x_{1}x_{2}\cdots
x_{n}}=\sum_{\hat{\mathcal{P}}}\int_{\ell_{2}\cdots\ell_{n}}\Psi^{A_{1}A_{2}\cdots
A_{n}}_{\ell_{2}\ell_{3}\cdots\ell_{n}|\hat{\mathcal{P}}}e^{-i(k_{\mathcal{P}}-\ell_{2\cdots
n)})\cdot x_{1}-\sum_{i=2}^{n}\ell_{i}\cdot x_{i}}\,,$ (60)

where $\ell_{12\cdots n}=\ell_{1}+\ell_{2}+\cdots+\ell_{n}$. For the two-loop
computations, we need up to second-descendant fields

$\displaystyle\psi^{AB}_{xy}$
$\displaystyle=\sum_{\hat{\mathcal{P}}}\int_{\ell}\Psi_{\ell|\hat{\mathcal{P}}}^{AB}e^{-i(k_{\mathcal{P}}-p)\cdot
x-ip\cdot y}\,,$ (61) $\displaystyle\psi^{ABC}_{xyz}$
$\displaystyle=\sum_{\hat{\mathcal{P}}}\int_{\ell_{1}\ell_{2}}\Psi_{\ell_{1}\ell_{2}|\hat{\mathcal{P}}}^{ABC}e^{-i(k_{\mathcal{P}}-\ell_{12})\cdot
x-i\ell_{1}\cdot y-i\ell_{2}\cdot z}\,.$

As we have defined $\hbar$-expansion of the classical fields in (35), we may
expand the off-shell currents as

$\displaystyle\Phi^{I}_{\mathcal{P}}$
$\displaystyle=\sum_{n=0}^{\infty}\int_{\ell_{1}\ell_{2}\cdots\ell_{n}}\Phi^{{\scriptscriptstyle(n)}}_{\ell_{1}\cdots\ell_{n}}\big{|}_{\mathcal{P}}^{I}\bigg{(}\frac{\hbar}{i}\bigg{)}^{n}\,,$
(62) $\displaystyle J^{M}_{\mathcal{P}}$
$\displaystyle=\sum_{n=0}^{\infty}\int_{\ell_{1}\ell_{2}\cdots\ell_{n}}J^{{\scriptscriptstyle(n)}}_{\ell_{1}\cdots\ell_{n}}\big{|}_{\mathcal{P}}^{M}\bigg{(}\frac{\hbar}{i}\bigg{)}^{n}\,,$
$\displaystyle\Psi^{AB_{1}\cdots B_{m}}_{\ell_{1}\ell_{2}\cdots\ell_{m}}$
$\displaystyle=\sum_{n=0}^{\infty}\int_{\ell^{\prime}_{1}\ell^{\prime}_{2}\cdots\ell^{\prime}_{n}}\Psi^{{\scriptscriptstyle(n)}}_{\ell^{\prime}_{1}\cdots\ell^{\prime}_{n}}\big{|}^{A,B_{1},\cdots,B_{m}}_{\ell_{1}\cdots\ell_{m}|\hat{\mathcal{P}}}\bigg{(}\frac{\hbar}{i}\bigg{)}^{n}\,,$

where $\ell_{1},\ell_{2}\cdots$ are the loop momenta. This $\hbar$-expansion
is crucial in constructing a self-consistent recursion relation by truncating
unnecessary higher descendant currents.

###  3.2 Bracket and Descend operator

Let us consider a classical field $A_{x}$ that is expressed as the product of
other classical fields $B^{1}_{x},B^{2}_{x},\cdots,B^{n}_{x}$,

$A_{x}=B^{1}_{x}B^{2}_{x}\cdots B^{n}_{x}\,,$ (63)

and assume that their perturbiner expansions are given by

$A_{x}=\sum_{\mathcal{P}}\mathcal{A}_{\mathcal{P}}e^{-ik_{\mathcal{P}}\cdot
x}\,,\qquad
B^{i}_{x}=\sum_{\mathcal{P}}\mathcal{B}^{i}_{\mathcal{P}}e^{-ik_{\mathcal{P}}\cdot
x}\,,$ (64)

where $\mathcal{A}_{\mathcal{P}}$ and $\mathcal{B}^{i}_{\mathcal{P}}$ denote
the off-shell currents. By substituting these perturbiner expansion into (63),
we can express $\mathcal{A}_{\mathcal{P}}$ in terms of $\mathcal{B}^{i}$

$\mathcal{A}_{\mathcal{P}}=\sum_{\mathcal{P}=\mathcal{Q}_{1}\cup\mathcal{Q}_{2}\cdots\cup\mathcal{Q}_{n}}\mathcal{B}^{1}_{\mathcal{Q}_{1}}\mathcal{B}^{2}_{\mathcal{Q}_{2}}\cdots\mathcal{B}^{n}_{\mathcal{Q}_{n}}\,.$
(65)

A typical form of off-shell recursions involves cumbersome products and
summations over the off-shell currents. For concise expressions, we introduce
a multilinear bracket denoted as $\lceil\bullet\rfloor_{\mathcal{P}}$, which
effectively replaces the summations over products of the currents,

$\displaystyle\big{\lceil}\mathcal{A},\mathcal{B}\big{\rfloor}_{\mathcal{P}}=\sum_{\mathcal{P}=\mathcal{Q}\cup\mathcal{R}}\mathcal{A}_{\mathcal{Q}}\mathcal{B}_{\mathcal{R}}\,,\qquad\big{\lceil}A,B,C\big{\rfloor}_{\mathcal{P}}=\sum_{\mathcal{P}=\mathcal{Q}\cup\mathcal{R}\cup\mathcal{S}}A_{\mathcal{Q}}B_{\mathcal{R}}C_{\mathcal{S}}\,.$
(66)

If descendant currents are positioned in one of the entries of the bracket,
the summation over the words is replaced with the hatted words, including the
zero modes. By definition, the bracket is orderless,
$\big{\lceil}\mathcal{A},\mathcal{B}\big{\rfloor}_{\mathcal{P}}=\big{\lceil}\mathcal{B},\mathcal{A}\big{\rfloor}_{\mathcal{P}}$,
and satisfies the multilinear property,
$\big{\lceil}a\mathcal{A}+b\mathcal{B},\mathcal{C}\big{\rfloor}_{\mathcal{P}}=a\big{\lceil}\mathcal{A},\mathcal{C}\big{\rfloor}_{\mathcal{P}}+b\big{\lceil}\mathcal{B},\mathcal{C}\big{\rfloor}_{\mathcal{P}}$.
Using this bracket, (65) can be rewritten by

$\mathcal{A}_{\mathcal{P}}=\big{\lceil}\mathcal{B}^{1},\mathcal{B}^{2},\cdots\mathcal{B}^{n}\big{\rfloor}_{\mathcal{P}}\,.$
(67)

One may insert momenta at any slot with an off-shell current, and the momenta
inherits the word of the current in the same slot, for instance,

$\displaystyle\big{\lceil}k^{\rho}J^{M},k^{\sigma}J^{N}\big{\rfloor}_{\mathcal{P}}$
$\displaystyle=\sum_{\mathcal{P}=\mathcal{Q}\cup\mathcal{R}}k^{\rho}_{\mathcal{Q}}J^{M}_{\mathcal{Q}}k^{\sigma}_{\mathcal{R}}J^{N}_{\mathcal{R}}\,,$
(68)
$\displaystyle\big{\lceil}J^{M},k^{\rho}k^{\sigma}J^{N}\big{\rfloor}_{\mathcal{P}}$
$\displaystyle=\sum_{\mathcal{P}=\mathcal{Q}\cup\mathcal{R}}J^{M}_{\mathcal{Q}}k^{\rho}_{\mathcal{R}}k^{\sigma}_{\mathcal{R}}J^{N}_{\mathcal{R}}\,.$

We may also add a superscript labeling the $\hbar$-order to the bracket, which
is the total $\hbar$-order of the currents in the bracket

$\displaystyle\big{\lceil}\mathcal{A},\mathcal{B}\big{\rfloor}^{{\scriptscriptstyle(n)}}_{\mathcal{P}}$
$\displaystyle=\sum_{m=0}^{n}\sum_{\mathcal{P}=\mathcal{Q}\cup\mathcal{R}}\mathcal{A}^{{\scriptscriptstyle(n-m)}}_{\mathcal{Q}}\mathcal{B}^{{\scriptscriptstyle(m)}}_{\mathcal{R}}\,,$
(69)
$\displaystyle\big{\lceil}\mathcal{A},\mathcal{B},\mathcal{C}\big{\rfloor}^{{\scriptscriptstyle(n)}}_{\mathcal{P}}$
$\displaystyle=\sum_{m,p=0\atop n\geq
m+p}^{n}\sum_{\mathcal{P}=\mathcal{Q}\cup\mathcal{R}\cup\mathcal{S}}\mathcal{A}^{{\scriptscriptstyle(n-m-p)}}_{\mathcal{Q}}\mathcal{B}^{{\scriptscriptstyle(m)}}_{\mathcal{R}}\mathcal{C}^{{\scriptscriptstyle(p)}}_{\mathcal{S}}\,.$

We next define a differential operator $\Delta^{A}_{\mathcal{P}}$ and call it
the descend operator. It maps an off-shell current and descendant currents as
follows:

$\displaystyle\Delta^{A}_{p}\Phi^{I}_{\mathcal{P}}$
$\displaystyle=\Psi^{IA}_{p|\mathcal{P}}\,,$ (70)
$\displaystyle\Delta^{A}_{p}J^{M}_{\mathcal{P}}$
$\displaystyle=\Psi^{MA}_{p|\mathcal{P}}\,,$
$\displaystyle\Delta^{A}_{p}\Psi^{BC_{1}\cdots C_{n}}_{q_{1}\cdots
q_{n}|\hat{\mathcal{P}}}$ $\displaystyle=\Psi^{BC_{1}\cdots
C_{n}A}_{q_{1}\cdots q_{n}p|\hat{\mathcal{P}}}\,.$

We also denote the successive applications of the descend operator as

$\Delta^{A_{1},A_{2},\cdots,A_{n}}_{\ell_{1},\ell_{2},\cdots,\ell_{n}}\big{[}\bullet\big{]}=\Delta^{A_{n}}_{\ell_{n}}\bigg{[}\cdots\Delta^{A_{2}}_{\ell_{2}}\Big{[}\Delta^{A_{1}}_{\ell_{1}}\big{[}\bullet\big{]}\Big{]}\cdots\bigg{]}\,.$
(71)

Since $\Delta^{A}_{p}$ is a differential operator, it should satisfy the
Leibniz rule inside the bracket

$\displaystyle\Delta^{A}_{p}\big{\lceil}J^{M},J^{N}\big{\rfloor}_{\mathcal{P}}$
$\displaystyle=\big{\lceil}\Delta^{A}_{p}J^{M},J^{N}\big{\rfloor}_{\mathcal{P}}+\big{\lceil}J^{M},\Delta^{A}_{p}J^{N}\big{\rfloor}_{\mathcal{P}}\,,$
(72)
$\displaystyle=\big{\lceil}\Psi^{MA}_{p},J^{N}\big{\rfloor}_{\mathcal{P}}+\big{\lceil}J^{M},\Psi^{NA}_{p}\big{\rfloor}_{\mathcal{P}}\,,$
$\displaystyle=\sum_{\mathcal{P}=\hat{\mathcal{Q}}\cup\hat{\mathcal{R}}}\Big{(}\Psi^{MA}_{p|\hat{\mathcal{Q}}}J^{N}_{\hat{\mathcal{R}}}+J^{M}_{\hat{\mathcal{Q}}}\Psi^{NA}_{p|\hat{\mathcal{R}}}\Big{)}\,.$

Note that the summation in the last line should be for the hatted words. It
also acts on the external on-shell momenta $k^{\mu}_{\mathcal{P}}$, and one
can show that the external momenta are shifted

$\displaystyle\Delta^{A}_{p}\big{[}k^{\mu}_{\mathcal{P}}\big{]}=k^{\mu}_{\mathcal{P}}-p^{\mu}:=k^{\mu}_{\mathcal{P},p}\,.$
(73)

###  3.3 Quantum off-shell recursions up to two loops

We now derive the off-shell recursion relations by combining the DS equation
with the quantum perturbiner expansion up to two-loop order, following the
steps established in Lee:2022aiu . Initially, we substitute the perturbiner
expansions of the classical fields and their descendant fields into the DS
equations. After this substitution, the series is truncated to isolate
contributions up to the two-loop order. We then extract the coefficients for
an arbitrary rank $\mathcal{P}$ corresponding to the off-shell recursions
under investigation. The descendant currents that contribute to the two-loop
recursions are

$\Psi^{{\scriptscriptstyle(0)}}\big{|}^{AB}_{\ell|\hat{\mathcal{P}}}\,,\qquad\Psi^{{\scriptscriptstyle(1)}}\big{|}^{AB}_{\ell|\hat{\mathcal{P}}}\,,\qquad\Psi^{{\scriptscriptstyle(0)}}\big{|}^{ABC}_{\ell_{1}\ell_{2}|\hat{\mathcal{P}}}\,.$
(74)

It is straightforward to derive the off-shell recursions for these descendant
currents by acting on the $\Delta$ operator.

####  3.3.1 Scalar field sector

The DS equation and the perturbiner expansions for the scalar fields are
derived in (36) and (55) respectively. Substituting the expansion and
truncating the DS equation at rank $|\mathcal{P}|$, in our case
$\mathcal{P}=123$ and $|\mathcal{P}|=3$, we present the recursions for the
off-shell currents for the massive scalars up to two loops,

$\displaystyle\big{(}k_{\mathcal{P}}^{2}+m_{I}^{2}\big{)}\Phi^{{\scriptscriptstyle(n)}}\big{|}^{I}_{\mathcal{P}}$
$\displaystyle=\kappa\Omega^{{\scriptscriptstyle(n)}}_{1}\big{|}^{I}_{\mathcal{P}}-\kappa^{2}\mathcal{I}^{M,N}_{\phi}\Omega^{{\scriptscriptstyle(n)}}_{2}\big{|}^{MN}_{\mathcal{P}}-\kappa^{3}\mathcal{I}^{M,N,P}_{\phi}\Omega^{{\scriptscriptstyle(n)}}_{3}\big{|}^{MNP}_{\mathcal{P}}\,,\quad
n\leq 2\,,$ (75)

where $\Omega^{{\scriptscriptstyle(n)}}_{i}$ are denoted in (270).

For consistency, we should require that each
$\Omega^{{\scriptscriptstyle(n)}}_{i}$ satisfies

$\Omega^{{\scriptscriptstyle(n)}}_{i}=0\,,\qquad\text{for}~{}n<0\,.$ (76)

One can derive the off-shell recursions for the descendant currents for
$\Phi^{I}_{\mathcal{P}}$ by acting the $\Delta$-operators on (75) by using the
relations in (70) and (73).

####  3.3.2 Graviton sector

Let us now consider the recursions for graviton current $J^{M}_{\mathcal{P}}$
and its descendants. We derive the recursions by substituting the perturbiner
expansion into the DS equation for $h^{\mu\nu}$ in (37). Collecting the $n$-th
order terms in $\hbar$, where $n\leq 2$, and the coefficients of
$e^{ik_{\mathcal{P}}\cdot x}$ provides the off-shell recursion relation. The
structure of the graviton recursion relation is the same as the DS equation
(37). The recursion at the $n$-loop order is given by

$\displaystyle
k^{2}_{\mathcal{P}}J^{\scriptscriptstyle(n)}\big{|}_{\mathcal{P}}^{M}$
$\displaystyle=2P^{M,N}\Bigg{[}\sum_{i=1}^{2}\kappa^{i}\mathcal{G}^{{\scriptscriptstyle(n)}}_{i}\big{|}^{N}_{\mathcal{P}}+\sum_{i=1}^{3}\kappa^{i}\mathcal{T}^{{\scriptscriptstyle(n)}}_{i}\big{|}^{N}_{\mathcal{P}}\Bigg{]}\,.$
(77)

Here the curvature part
$\mathcal{G}^{{\scriptscriptstyle(n)}}_{i}\big{|}^{M}_{\mathcal{P}}$ are
defined as

$\displaystyle\mathcal{G}^{{\scriptscriptstyle(n)}}_{1}\big{|}^{N}_{\mathcal{P}}$
$\displaystyle=\tilde{\mathcal{I}}^{M,\rho N,\sigma
P}_{1}\Lambda^{\scriptscriptstyle(n)}_{1,1}\big{|}^{\rho N,\sigma
P}_{\mathcal{P}}+\tilde{\mathcal{I}}^{M,N,\rho\sigma,P}_{2}\Lambda^{\scriptscriptstyle(n)}_{1,2}\big{|}^{N,\rho\sigma,P}_{\mathcal{P}}\,$
(78)
$\displaystyle\mathcal{G}^{{\scriptscriptstyle(n)}}_{2}\big{|}^{N}_{\mathcal{P}}$
$\displaystyle=\tilde{\mathcal{I}}^{M,N,\rho P,\sigma
Q}_{1}\Lambda^{\scriptscriptstyle(n)}_{2,1}\big{|}^{N,\rho P,\sigma
Q}_{\mathcal{P}}+\tilde{\mathcal{I}}_{2}^{M,N,P,\rho\sigma,Q}\Lambda^{\scriptscriptstyle(n)}_{2,2}\big{|}^{N,P,\rho\sigma,Q}_{\mathcal{P}}\,,$

and the explicit form of the $\Lambda^{{\scriptscriptstyle(n)}}_{n,1}$ and
$\Lambda^{{\scriptscriptstyle(n)}}_{n,2}$ are listed in (271). The EM tensor
part $\mathcal{T}^{{\scriptscriptstyle(n)}}_{i}\big{|}^{N}_{\mathcal{P}}$ are
defined as

$\displaystyle\mathcal{T}^{{\scriptscriptstyle(n)}}_{1}\big{|}^{M}_{\mathcal{P}}$
$\displaystyle=\Theta_{1}^{\scriptscriptstyle(n)}\big{|}^{M}_{\mathcal{P}}\,,$
(79)
$\displaystyle\mathcal{T}^{{\scriptscriptstyle(n)}}_{2}\big{|}^{M}_{\mathcal{P}}$
$\displaystyle=\mathcal{I}_{\rm
EM}^{M,N}\Theta_{2}^{\scriptscriptstyle(n)}\big{|}^{N}_{\mathcal{P}}\,,$
$\displaystyle\mathcal{T}^{{\scriptscriptstyle(n)}}_{3}\big{|}^{M}_{\mathcal{P}}$
$\displaystyle=\mathcal{I}_{\rm
EM}^{M,N,P}\Theta_{3}^{\scriptscriptstyle(n)}\big{|}^{N,P}_{\mathcal{P}}\,,$

and the explicit form of $\Theta^{\scriptscriptstyle(n)}_{i}$ are listed in
(272).

Again, we require that $\Lambda^{{\scriptscriptstyle(n)}}$ and
$\Theta^{{\scriptscriptstyle(n)}}$ vanish for $n<0$ cases

$\Lambda_{i,j}^{{\scriptscriptstyle(n)}}=0\,,\qquad\Theta_{i}^{{\scriptscriptstyle(n)}}=0\,,\qquad\text{for}~{}n<0\,,$
(80)

and the associated descendent recursions can be derived by acting the
$\Delta$-operators successively to (77).

###  3.4 Initial conditions

We should impose an appropriate initial condition to address the solutions of
the off-shell recursions. For $2\to 2$ scattering amplitudes for scalar
particles, the external states consist of scalar fields only, while gravitons
are restricted to the internal lines. The initial condition for the recursions
is realized by the rank-1 off-shell currents for the scalars $\Phi^{I}_{i}$
and for graviton $J^{\mu\nu}_{i}$, as well as the descendant fields.

First, let us consider $\Phi^{I}_{i}$ and $J^{\mu\nu}_{i}$. By definition of
the off-shell currents and the dressed external source, (81) and (43), we may
derive the initial condition

$\displaystyle\Phi^{I}_{i}$
$\displaystyle=\big{\langle}0\big{|}\Phi^{I}(-k_{i})\Phi^{J}(k_{i})\big{|}0\big{\rangle}\tilde{\mathbf{K}}^{J,K}_{k_{i}}\epsilon^{K}_{i}=\tilde{\mathbf{D}}^{I,J}_{k_{i}}\tilde{\mathbf{K}}^{J,K}_{k_{i}}\epsilon^{K}_{i}=\epsilon^{I}_{i}\,,$
(81) $\displaystyle J^{\mu\nu}_{i}$ $\displaystyle=0\,.$

Note that this initial condition is for the all-loop order, and we should
distribute the result to each loop order. In fact, one can show that the
initial condition is only contributed from the rank-1 currents at the tree
level, and the higher loops vanish. The tree-level two-point function, which
is equivalent to the propagator $\tilde{D}^{I,J}_{k_{i}}$, and the inverse
propagator $\tilde{K}^{I,J}_{k_{i}}$ satisfy the inverse relation
$\tilde{D}^{I,J}_{k_{i}}\tilde{K}^{J,K}_{k_{i}}=\delta^{IK}$.

Consequently, the first equation of (81) can be expanded in $\hbar$ as

$\tilde{\mathbf{D}}^{I,J}_{k_{i}}\tilde{\mathbf{K}}^{J,K}_{k_{i}}=\delta^{IK}+\sum_{n=1}^{\infty}\sum_{m=1}^{n}\Big{(}\frac{\hbar}{i}\Big{)}^{n}\tilde{\mathbf{D}}^{{\scriptscriptstyle(n-m)}}\big{|}^{I,J}_{k_{i}}\tilde{\mathbf{K}}^{{\scriptscriptstyle(m)}}\big{|}^{J,K}_{k_{i}}=\delta^{IK}\,.$
(82)

As this equation holds for each order in $\hbar$, contributions from higher
loop orders vanish

$\Phi^{{\scriptscriptstyle(n)}}|^{I}_{i}=\sum_{m=1}^{n}\Big{(}\frac{\hbar}{i}\Big{)}^{n}\tilde{\mathbf{D}}^{{\scriptscriptstyle(n-m)}}\big{|}^{I,J}_{k_{i}}\tilde{\mathbf{K}}^{{\scriptscriptstyle(m)}}\big{|}^{J,K}_{k_{i}}\epsilon^{K}=0\,,\qquad\text{for}~{}n>0\,.$
(83)

Thus our initial condition for $\Phi^{I}_{\mathcal{P}}$ and
$J^{\mu\nu}_{\mathcal{P}}$ is

$\displaystyle\Phi^{{\scriptscriptstyle(0)}}\big{|}^{I}_{i}$
$\displaystyle=\epsilon^{I}_{i}\,,$ $\displaystyle\qquad
J^{{\scriptscriptstyle(0)}}\big{|}^{I}_{i}$ $\displaystyle=0\,,$ (84)
$\displaystyle\Phi^{{\scriptscriptstyle(n)}}\big{|}^{I}_{i}$
$\displaystyle=0\,,$ $\displaystyle\qquad
J^{{\scriptscriptstyle(n)}}\big{|}^{I}_{i}$
$\displaystyle=0\,,\qquad\text{for}~{}n>0\,.$

The DS equations with the bare external sources are reduced to

$\displaystyle\sum_{i}k_{i}^{2}\big{(}J^{{\scriptscriptstyle(0)}}|^{\mu\nu}_{i}e^{-ik_{i}\cdot
x}\big{)}$
$\displaystyle=4P^{\mu\nu,\rho\sigma}j^{\rho\sigma}_{x}=\sum_{i}k^{2}\epsilon^{\mu\nu}_{i}e^{ik_{i}\cdot
x}\,,$ (85)
$\displaystyle\sum_{i}\big{(}k_{i}^{2}+m^{2}_{I}\big{)}\big{(}\Phi^{I}_{i}e^{-ik_{i}\cdot
x}\big{)}$
$\displaystyle=\rho^{I}_{x}=\big{(}k_{i}^{2}+m^{2}_{I}\big{)}\epsilon^{I}_{i}e^{-ik_{i}\cdot
x}\,,$

and the previous initial conditions satisfy them.

Finally, let us consider the initial condition for the descendant currents
case. As before, we derive the initial condition from their definition. One
can show that it is enough to determine the case of the first descendant
currents because the recursion relations generate the initial conditions for
the other descendant currents. Unlike the $J^{\mu\nu}_{\mathcal{P}}$ and
$\Phi^{I}_{\mathcal{P}}$, the descendant currents may have the zero-mode
contribution

$\displaystyle\Psi^{AB}_{\ell,\emptyset}=\langle
0|\tilde{\mathbf{h}}^{A}(-\ell)\tilde{\mathbf{h}}^{B}(\ell)|0\rangle_{c}=\tilde{\mathbf{D}}^{AB}_{\ell}\,.$
(86)

However, the higher loop currents
$\Psi^{{\scriptscriptstyle(n>0)}}\big{|}^{AB}_{\ell,\emptyset}$ contribute to
the self-energy or the vacuum polarizations, which are classified as quantum
contributions according to the prescription discussed in Section 4.2. Thus
only the tree-level current
$\Psi^{{\scriptscriptstyle(0)}}\big{|}^{AB}_{\ell,\emptyset}$ contributes to
the classical amplitudes. Thus, our initial condition is

$\displaystyle\Psi^{{\scriptscriptstyle(0)}}\big{|}^{IJ}_{\ell,\emptyset}$
$\displaystyle=\frac{\delta^{IJ}}{\ell^{2}+m^{2}_{I}}\,,$
$\displaystyle\qquad\Psi^{{\scriptscriptstyle(n)}}\big{|}^{IJ}_{\ell,\emptyset}$
$\displaystyle=0\,,\quad\text{for}~{}n>0\,,$ (87)
$\displaystyle\Psi^{{\scriptscriptstyle(0)}}\big{|}^{MN}_{\ell,\emptyset}$
$\displaystyle=\frac{2P^{M,N}}{\ell^{2}}\,,$
$\displaystyle\qquad\Psi^{{\scriptscriptstyle(n)}}\big{|}^{MN}_{\ell,\emptyset}$
$\displaystyle=0\,,\quad\text{for}~{}n>0\,,$

and these are propagators and represented diagrammatically in Figure 1.

$\Psi^{\mu\nu,\rho\sigma}_{p,\emptyset}$$=$$-p$$\mu\nu$$p$$\rho\sigma$Tree$=$$p$$\mu\nu$$\rho\sigma$
Figure 1: the wiggly lines represent off-shell gravitons

One can show that
$\Psi^{{\scriptscriptstyle(0)}}\big{|}^{AB}_{\ell,\emptyset}$ satisfies the DS
equation as we expect:

$\displaystyle\int_{\ell}\left(\square_{x}-m_{I}^{2}\right)\big{(}\Psi^{{\scriptscriptstyle(0)}}\big{|}^{IJ}_{\ell|\emptyset}e^{i\ell\cdot(x-y)}\big{)}$
$\displaystyle=-\int_{\ell}\delta^{IJ}e^{i\ell\cdot(x-y)}\,,$ (88)
$\displaystyle\int_{\ell}\square_{x}\big{(}\Psi^{{\scriptscriptstyle(0)}}\big{|}^{MN}_{\ell|\emptyset}e^{i\ell\cdot(x-y)}\big{)}$
$\displaystyle=-\int_{\ell}2P^{MN}e^{i\ell\cdot(x-y)}\,.$

Interestingly, all the other descendant currents are determined by the
$\Phi^{I}_{i}$, $J^{\mu\nu}_{i}$ and $\Psi^{AB}_{\ell,\emptyset}$ by the
recursion relations. Consequently, the total initial condition is given by
(84) and (87).

##  4 Classifying Classical Diagrams

So far, we have constructed the quantum off-shell recursions for the
scattering amplitude of scalar particles up to the two-loop order. We may
solve the recursions using symbolic computational tools, such as Mathematica,
and the solutions of the off-shell recursions serve as loop integrands. The
loop integrals include not only quantum but also classical contributions.
Since the quantum effects are negligible in the binary black hole systems, we
will focus on establishing a method for identifying the classical diagrams
pertinent to the off-shell recursions.

Notably, the off-shell recursion method offers an algebraic approach that
contrasts with the conventional unitarity cut method – it does not require any
explicit diagram topologies for solving the recursions. Consequently, we
intend to develop a novel prescription to identify and extract the classical
diagrams by employing power counting rules, which do not rely on explicit
diagrammatic representation.

###  4.1 Kinematics

$k_{2},m_{2}$$k_{3},m_{2}$$k_{1},m_{1}$$k_{4},m_{1}$ Figure 2: The kinematics
of binary BH systems

Let us consider the classical $2\to 2$ scattering amplitude $\mathcal{M}_{\rm
cl}[k_{1},k_{2},k_{3},k_{4}]$ of two massive scalar particles, where $k_{1}$,
$k_{2}$, $k_{3}$ and $k_{4}$ are external momenta for the scalar particles as
depicted in Figure 2. From the relation between the off-shell current and the
amplitude (56), the amplitude is associated with the rank-3 scalar off-shell
current $\Phi^{1}_{123}$ depicted in Figure 3. Note that $\Phi^{1}_{234}$ or
$\Phi^{2}_{124}$ etc would give the same result, but here we will focus on
$\Phi^{1}_{123}$ only.

The external momenta satisfy the momentum conservation,
$k_{1}+k_{2}+k_{3}+k_{4}=0$, and the mass-shell condition,
$k_{1}^{2}+m_{1}^{2}=0$, $k_{2}^{2}+m_{2}^{2}=0$, $k_{3}^{2}+m_{2}^{2}=0$,
$k_{4}^{2}+m_{1}^{2}=0$. For the soft expansion, which we will discuss in the
next section, it is useful to parametrize the external momenta Parra-
Martinez:2020dzs  as follows:

$k^{\mu}_{1}=-\bar{k}^{\mu}_{1}+\frac{q^{\mu}}{2}\,,\qquad
k^{\mu}_{2}=-\bar{k}^{\mu}_{2}-\frac{q^{\mu}}{2}\,,\qquad
k^{\mu}_{3}=\bar{k}^{\mu}_{2}-\frac{q^{\mu}}{2}\,,\qquad
k^{\mu}_{4}=\bar{k}^{\mu}_{1}+\frac{q^{\mu}}{2}\,,$ (89)

where $q^{\mu}$ is the momentum transfer between the scalar particles. The
mass-shell conditions require $\bar{k}^{\mu}_{1}$ and $\bar{k}^{\mu}_{2}$ are
orthogonal to $q^{\mu}$

$k^{2}_{1}-k^{2}_{4}=-2\bar{k}_{1}\cdot q=0\,,\qquad
k^{2}_{2}-k^{2}_{3}=2\bar{k}_{2}\cdot q=0\,.$ (90)

We also introduce the dimensionless momenta $u^{\mu}_{1}$ and $u^{\mu}_{2}$

$u^{\mu}_{1}=\frac{\bar{k}^{\mu}_{1}}{\bar{m}_{1}}\,,\qquad
u^{\mu}_{2}=\frac{\bar{k}^{\mu}_{2}}{\bar{m}_{2}}\,,$ (91)

where

$\bar{m}^{2}_{1}=-\bar{k}^{2}_{1}=m^{2}_{1}+\frac{q^{2}}{4}\,,\qquad\bar{m}^{2}_{2}=-\bar{k}^{2}_{2}=m^{2}_{2}+\frac{q^{2}}{4}\,.$
(92)

It is useful to introduce parameters $y$ and $\sigma$ as

$y=-u_{1}\cdot u_{2}\,,\qquad\sigma=-\frac{k_{1}\cdot k_{2}}{m_{1}m_{2}}\,,$
(93)

and these are related to

$y=\sigma-\frac{2m_{1}m_{2}+\sigma(m_{1}^{2}+m_{2}^{2})}{8m_{1}^{2}m_{2}^{2}}q^{2}+\cdots\,.$
(94)

We further introduce a parameter $x$ to rationalize $\sqrt{y^{2}-1}$ which
appears in the middle of computations

$x=y-\sqrt{y^{2}-1}\,,\qquad\sqrt{y^{2}-1}=\frac{1-x^{2}}{x}\,,\qquad y\geq
1\,,\qquad 0<x\leq 1\,.$ (95)

$=$$\Phi^{1}_{123}$ $k_{1}$$k_{2}$$k_{3}$ Figure 3: Diagrammatical
representation of the currents associated with the tree and one-loop scalar
four-point amplitude. The double lines represent the off-shell lines

###  4.2 Classical loop integrands

Solving the quantum off-shell recursions, we can derive the off-shell current
$\Phi^{1}_{123}$ associated with the $2\to 2$ scattering amplitude for the
scalar particles at a given loop order. Interestingly, even at the loop
levels, the amplitude includes both classical and quantum contributions. For
describing the classical conservative dynamics, we only focus on the classical
part.

For a consistent field theory limit, we assume that the black holes are well
separated to ignore the internal structure of the black holes, and the Compton
length scale must be much smaller than the impact parameter $|b|$. The Fourier
dual of $|b|$ is given by the momentum transfer $|q|$ between two black holes,
$|b|\sim\frac{1}{|q|}$, and this indicates $|q|$ is much smaller than the mass
and the momenta of the black holes

$m_{1},m_{2},|k_{i}|\gg|q|\,.$ (96)

This implies that we may characterize the amplitude according to the expansion
in small $|q|$. We further consider the conservative dynamics by assuming the
small relative velocity limit $\vec{v}\ll 1$, or a scale of the components of
$q^{\mu}=(q^{0},\vec{q})$

$q^{0}\gg\vec{q}\,.$ (97)

Thus we consider two limits – near-static classical limit, $|q|\ll 1$ and
$\vec{v}\ll 1$.

To isolate the classical contributions from the total scattering amplitude, it
is necessary to identify the classical loop integrands contributing to the
classical amplitudes. A simple and efficient prescription was introduced in
Bern:2019crd . The established criteria for a classical loop integrands are:

  * •

Every independent loop must contain at least one matter line

  * •

Both ends of a graviton line cannot be connected to the same matter line

  * •

Truncate the higher power of loop momenta and momentum transfer

The first two are rules on the topology of the diagrams, and the diagrams
satisfying the criteria are called classical diagrams. The last one is a
condition on the numerators of loop integrands, which is realized by the soft
expansion in Section 5. This prescription is particularly useful in the
context of the unitarity cut method because cut topologies should be imposed
before calculation thereby one can focus on the classical diagrams only.

In contrast, the quantum off-shell recursion formalism generates the loop
integrands without concerning the diagram topologies by solving the recursions
algebraically. However, it is challenging to figure out the diagram topologies
from the total loop integrands. Thus we need an alternative strategy for
identifying classical diagrams that do not rely on diagram topologies.

To address this issue, we propose a novel prescription that employs power
counting rules to separate classical diagrams. To this end, we introduce
additional fictitious coupling constants and labeling for each propagator. The
only coupling constant in the field theory is the Newton constant, $G_{N}$ or
$\kappa$, and the loop integrands at $L$-loop order have the same couplings
dependence, $\kappa^{2(L+1)}$. For a finer separation, we introduce new formal
coupling constants, $\kappa_{1}$, $\kappa_{2}$ and $\kappa_{g}$, obtained by
dividing $\kappa$ according to the types of interaction vertice as depicted in
Figure 4. These are useful bookkeeping devices that convey more information
than $\kappa$.

matter 1:$(\kappa_{1})^{1}$

$(\kappa_{1})^{2}$

$(\kappa_{1})^{3}$

matter 2:$(\kappa_{2})^{1}$

$(\kappa_{2})^{2}$

$(\kappa_{2})^{3}$

graviton:$(\kappa_{g})^{1}$

$(\kappa_{g})^{2}$

Figure 4: We replace the Newton constant $\kappa$ to new formal coupling
constants $\kappa_{1}$, $\kappa_{2}$ and $\kappa_{g}$ according to the
interaction types.

Next, we introduce labels for each propagator, $n_{1}$, $n_{2}$ and $n_{g}$,
for the $\phi_{1}$, $\phi_{2}$ and the graviton respectively

graviton : $n_{g}$matter 1 : $n_{1}$matter 2 : $n_{2}$ (98)

Collecting the couplings and labelings, we can define an identifier of
diagrams as follows:

$(\kappa_{g})^{N_{1}}(\kappa_{1})^{N_{2}}(\kappa_{2})^{N_{3}}(n_{g})^{N_{4}}(n_{1})^{N_{5}}(n_{2})^{N_{6}}\,.$
(99)

We may assign a set of the six numbers $(N_{1},\cdots N_{6})$ for each
diagram. These numbers can be facilitated to define conditions on the
identifiers for filtering the classical diagrams efficiently. However, this
assignment is not an one-to-one correspondence between identifiers and
diagrams. Multiple diagrams can correspond to a set of exponents in (99); even
quantum diagrams can satisfy the conditions for classical diagrams.
Nevertheless, such quantum diagrams can be effectively removed because of
their distinctive propagator structure compared to classical diagrams. This is
sufficient for our purpose.

There is a subtle issue in assigning the identifier. Let us consider a diagram
having graviton propagators without a loop momentum, which involves
$\frac{1}{q^{2}}$ factor. Since the factor does not affect the loop integrals,
we may consider them as diagrams without such a graviton propagator. While
these diagrams may exhibit distinct power-counting behaviors, we do not treat
them as independent diagrams as in the following examples,

$=\frac{1}{q^{2}}$ $=\frac{1}{q^{2}}$ (100)

This feature simplifies the classification scheme discussed in the following
subsection.

In practice, we do not have to impose those numbers by hand. Instead, we can
assign the numbers at the level of recursions. It is achieved simply by
replacing $\kappa$ to $\\{\kappa_{1},\kappa_{1},\kappa_{g}\\}$ and scaling
each propagator in the recursions such as

$\frac{1}{(\ell\pm k_{\mathcal{P}})^{2}}\to\frac{n_{g}}{(\ell\pm
k^{2}_{\mathcal{P}})^{2}}\,,\qquad\frac{1}{(\ell\pm
k_{\mathcal{P}})^{2}+m^{2}_{1,2}}\to\frac{n_{1,2}}{(\ell\pm
k_{\mathcal{P}})^{2}+m^{2}_{1,2}}\,.$ (101)

For instance, we may rewrite the recursion of $\Phi^{I}$ in (75) as

$\displaystyle\Phi^{{\scriptscriptstyle(2)}}_{\ell_{1}\ell_{2}}\big{|}^{I}_{\mathcal{P}}=\frac{n_{I}}{k_{\mathcal{P}}+m^{2}_{I}}\Big{(}\kappa_{I}\Omega^{{\scriptscriptstyle(2)}}_{\ell_{1}\ell_{2}}\big{|}^{I}_{\mathcal{P}}-\kappa_{I}^{2}\mathcal{I}^{M,N}_{\varphi,I}\Omega^{{\scriptscriptstyle(2)}}_{\ell_{1}\ell_{2}}\big{|}^{I,M,N}_{\mathcal{P}}-\kappa_{I}^{3}\mathcal{I}^{M,N,P}_{\varphi,I}\Omega^{{\scriptscriptstyle(2)}}_{\ell_{1}\ell_{2}}\big{|}^{I,M,N,P}_{\mathcal{P}}\Big{)}\,.$
(102)

If we solve the modified recursions, each term in the off-shell currents
contains the polynomial of the couplings and the labeling of the propagators
automatically.

###  4.3 Classifying the classical diagrams

We now introduce a new classification scheme for efficient evaluations of loop
integrals, especially optimizing the integration-by-parts (IBP) reduction
procedure. Before discussing IBP reduction in the next section, we need to
introduce integral family for a diagram.

Typically, one may associate an integral family for each diagram. An integral
family $F^{\Gamma_{l}}_{i_{1},i_{2},\cdots,i_{n}}$ at $l$ loops associated
with a diagram $\Gamma_{l}$ is represented by

$F^{\Gamma_{l}}_{i_{1},i_{2},\cdots,i_{n}}=\int\frac{\mathrm{d}^{D}\ell_{1}}{(2\pi)^{D}}\cdots\frac{\mathrm{d}^{D}\ell_{l}}{(2\pi)^{D}}\frac{1}{\big{(}\rho_{1}\big{)}^{i_{1}}\big{(}\rho_{2}\big{)}^{i_{2}}\cdots\big{(}\rho_{n}\big{)}^{i_{n}}}\,,\qquad
i_{1},\cdots,i_{n}\in\mathbb{Z}\,,$ (103)

where $\rho_{i}$ represent the propagators that constitute the diagram
$\Gamma_{l}$. Since the exponents can be zero or negative integers, the above
is the minimum form – one may add several propagators not included in the loop
integrand associated with the diagram $\Gamma_{l}$. Note that the assignment
of integral families for each diagram is computationally redundant.

Instead, our method offers a more efficient alternative: grouping diagrams
characterized by the integral family. After clustering the diagrams, the IBP
reduction can be executed on the groups collectively rather than individually.
This collective approach significantly enhances the computational efficiency.

####  4.3.1 Classes and Sectors

Our classification scheme occurs through a two-step process. The first step is
grouping the classical diagrams according to the power counting of the new
couplings $\kappa_{1},\kappa_{2}$ and $\kappa_{g}$. To facilitate this, we
introduce a process called pinching, which removes either an internal matter
or a graviton line (propagator) from a diagram. While this process reduces the
total number of internal lines, it leaves the powers of the coupling constants
$\kappa_{1},\kappa_{2}$ and $\kappa_{g}$ unaltered. We may illustrate the
pinching process as follows:

pinching (104) pinching

Thus, any integrands related through the pinching share the same powers of
coupling constants. This invariance enables us to define a specific set of
diagrams, which we shall call this collection as a class. All the integrands
constituting a given class carry the same power-counting of the coupling
constants regardless of the number of internal lines. Thus, the powers of the
couplings can be used for markers of classes:

$\text{class}\iff\big{(}\kappa_{g}\big{)}^{N_{1}}\big{(}\kappa_{1}\big{)}^{N_{2}}\big{(}\kappa_{2}\big{)}^{N_{3}}\,.$
(105)

The next step is classifying according to the structure of propagators, in
other words, the diagram topologies. The power-counting criterion alone is
insufficient for distinguishing the difference among all diagram topologies.
This implies that a class allows a further decomposition into multiple
subclasses according to the available integral families, which depend on the
structure of propagators. We shall call such subclasses sector. Obviously, the
choice of the subclasses is not unique. We want a systematic construction that
minimizes the number of subclasses.

Let us consider specific diagrams with maximal internal lines, which maximize
the sum of $n_{1}+n_{2}+n_{g}$ within a class and refer to these as primary
diagrams. In organizing sectors within a class, we use these primary diagrams
as a reference point. More precisely, a sector comprises a collection of
diagrams directly related to a designated primary diagram through pinching
while keeping the structure of the interaction vertices. Generally, the
primary diagram is defined uniquely for a given class, and we may use the
primary diagram to label the class.

A salient feature of the sector is that every integrand within a given sector
belongs to the same integral family. This consistency arises from the property
of pinching: it removes an internal propagator by setting the power of the
propagator to zero, yet the resultant integrand remains in the original
integral family. Consequently, each primary diagram defines a sector, and each
sector is associated with a specific integral family. In Figure 5, we
summarized the sector.

ClassSector 2Sector 1Sector 3Integral Family 2Integral Family 1Integral Family
3 Figure 5: A class is decomposed by each sector, generated from a primary
diagram by pinching. An integral family is associated with each sector.

One important remark is that constructing a sector is not unique. For
instance, consider the process of obtaining a four-point vertex through
pinching as below,

pinching (106)

This implies that a pinched diagram can have multiple origins, and this
feature creates ambiguity in defining a sector. However, this ambiguity does
not affect the final results. As such, the choice of a sector remains flexible
and can be adapted depending on our purpose.

####  4.3.2 Crossing and mirror duals

Let us consider a set of all possible diagrams associated with an amplitude.
For a given diagram, one can always identify its crossing and mirror dual
counterparts within the set. Combining these diagrams, we can define a
quartet. The loop integrands associated with the elements of a quartet are
related to each other in a specific way. This structure enables us to reduce
the number of loop integrals to be performed.

IY

$\overline{\rm IY}$

YI

$\rm IY$

Figure 6: A quartet including the IY diagram, which includes its crossed, the
horizontal mirror dual and the vertical mirror dual diagrams are shown in
order.

To begin with, we explore the concept of crossing. An integrand for a crossed
diagram is generated by swapping the momenta $k_{2}$ and $k_{3}$ for a given
diagram. This exchange is represented as follows:

$\displaystyle\text{\bf crossing}:~{}k_{2}\leftrightarrow
k_{3}\Longrightarrow{}$ $\displaystyle u_{1}\to u_{1}\,,~{}u_{2}\to-
u_{2}\,,~{}y\to-y\,,~{}q^{\mu}\to q^{\mu}\,,$ (107)
$\displaystyle\sigma\to-\sigma+\frac{q^{2}}{2m_{1}m_{2}}\,,~{}\sqrt{y^{2}-1}\to-\sqrt{y^{2}-1}\,,\quad
x\to-x\,.$

We denote the crossed diagram as $\overline{\Gamma}$ for a diagram $\Gamma$.
We illustrate an example in Figure 6. The second diagram ($\overline{\rm IY}$)
is the crossed diagram of the first diagram (IY). Note that the crossing map
modifies the pole structure of loop integrals at the potential region, and we
cannot be analytically continued from $y>1\to y<-1$ Parra-Martinez:2020dzs .
Thus, this relation should be applied at the integrand level rather than the
loop integral values.

Next, let us consider the mirror dual diagrams generated by reflecting the
spatial orientation of the original diagram. We distinguish between two types
of mirror duals, namely horizontal and vertical, formally defined as:

$\displaystyle\text{\bf horizontal mirror dual}:$
$\displaystyle~{}k_{1}\leftrightarrow
k_{4}~{}\text{and}~{}k_{2}\leftrightarrow k_{3}\Longrightarrow y\to
y\,,~{}q^{\mu}\to q^{\mu}\,,$ (108) $\displaystyle\text{\bf vertical mirror
dual}:$ $\displaystyle~{}k_{1}\leftrightarrow
k_{2}~{}\text{and}~{}k_{3}\leftrightarrow k_{4}\Longrightarrow y\to
y\,,~{}q^{\mu}\to-q^{\mu}\,,~{}m_{1}\leftrightarrow m_{2}\,.$

An example of the horizontal and vertical mirror duals is illustrated in
Figure 6. Unlike the crossing, these mirror duals preserve the sign of $y$,
circumventing any concerns related to the pole structure. The mirror dual
relations can be directly applied to the integral values, and we do not have
to evaluate loop integrals for mirror dual diagrams.

The total loop integrands, the solution of the recursions at a specific loop
order, contain all possible diagram topologies. Within the set of total
integrands, one may find the crossing and mirror dual pairs for a given
diagram, including the self-dual diagrams. According to our classification
scheme, we may extend the notion of crossing and mirror duals to the level of
sectors. This structure significantly enhances computational efficiency
because the symmetry properties inherent in the set of diagrams render it
unnecessary to evaluate the loop integrals corresponding to each mirror dual
independently. Furthermore, such relations among the sectors serve as a robust
consistency check for the loop integrands.

####  4.3.3 One-loop example

box (II)

crossed box ($\overline{\rm II}$)

triangle ($\triangle$)

inverted triangle (

$\triangle$

)

Figure 7: One loop diagrams

To elucidate the definitions previously outlined, such as the class, sector
and integral family, let us examine the classical $2\to 2$ amplitude at one
loop. It is well-known that the classical one-loop contribution arises from
both box and triangle diagrams Bern:2019crd ; Cristofoli:2020uzm . There exist
two kinds of box diagrams: the box (II) diagram and the crossed box
($\overline{\rm II}$), which are illustrated in figure 7. These have the same
power counting,

$\text{II, $\overline{\rm
II}$}:~{}(\kappa_{1})^{2}(\kappa_{2})^{2}(\kappa_{g})^{0}(n_{1})^{1}(n_{2})^{1}(n_{g})^{2}\,.$
(109)

Similarly, there are two kinds of triangle diagrams: the triangle
($\triangle$) and the inverted triangle (  $\triangle$ ). Their power
countings are given by222Note that we omitted the diagrams containing a
graviton propagator without loop momenta, $\frac{1}{q^{2}}$ as in (100).

$\displaystyle\triangle:{}$
$\displaystyle(\kappa_{1})^{2}(\kappa_{2})^{2}(\kappa_{g})^{0}(n_{1})^{0}(n_{2})^{1}(n_{g})^{2}\,,$
(110) $\displaystyle\bigtriangledown:{}$
$\displaystyle(\kappa_{1})^{2}(\kappa_{2})^{2}(\kappa_{g})^{0}(n_{1})^{1}(n_{2})^{0}(n_{g})^{2}\,.$

One can check that the box and the triangle diagrams have the same power-
counting of the couplings, $\kappa^{2}_{1}\kappa^{2}_{2}\kappa^{0}_{g}$.
According to the relation between a class and power-counting of the couplings
(105), all the classical diagrams at one-loop form a single class, which we
shall refer to $\Box$-class. Furthermore, the triangle diagrams can be
generated by pinching a matter line from the box diagrams.

The primary diagrams in $\Box$-class are the II and $\overline{\rm II}$
diagrams. These have four internal lines, whereas the triangle diagrams have
three lines. Each primary diagram defines a sector: the II-sector and the
$\overline{\rm II}$-sector. Thus, the entire one-loop integrands are
decomposed into the two sectors. As described in (106), there is no canonical
choice for defining a sector. We should determine how to allocate the triangle
diagrams into the sectors. Our convention is to equally distribute the
triangles into each sector as follows:

$\displaystyle\text{II-sector}:$
$\displaystyle~{}\mbox{II}+\frac{1}{2}\big{(}\triangle+\scalebox{-1.0}[1.0]{\rotatebox[origin={c}]{180.0}{$\triangle$}}\big{)}\,,$
(111) $\displaystyle\text{$\overline{\rm II}$-sector}:$
$\displaystyle~{}\mbox{X}+\frac{1}{2}\big{(}\triangle+\scalebox{-1.0}[1.0]{\rotatebox[origin={c}]{180.0}{$\triangle$}}\big{)}\,.$

Since the triangle diagrams are invariant under the crossing relation (107),
this choice ensures a manifest invariance under the crossing at the level of
integrands. We may exploit this property as a consistency check of the
solutions of the off-shell recursion.

There are two integral families associated with the II and the $\overline{\rm
II}$ sectors. First, the integral family for the II sector is defined by

$\displaystyle F^{\rm II}_{i_{1},i_{2},i_{3},i_{4}}$
$\displaystyle=\int_{\ell}\frac{1}{\big{(}\rho^{\rm
II,1}_{g}\big{)}^{i_{1}}\big{(}\rho^{\rm
II,2}_{g}\big{)}^{i_{2}}\big{(}\rho^{\rm
II,3}_{1}\big{)}^{i_{3}}\big{(}\rho^{\rm II,4}_{2}\big{)}^{i_{4}}}\,,\qquad
i_{1},\cdots i_{4}\in\mathbb{Z}\,,$ (112)

where the loop integral is defined by Smirnov:2012gma ; Parra-Martinez:2020dzs

$\int_{\ell}\equiv\int\frac{d^{D}\ell\,e^{\gamma_{E}\epsilon}}{i\pi^{D/2}}\,,$
(113)

which is different from the usual convention of the Fourier transformation,
$\int\frac{\mathrm{d}^{4}\ell_{i}}{(2\pi)^{4}}$. We should multiply the
following factor for each loop integral to restore the original convention and
compute the scattering amplitude,

$\frac{i}{\big{(}4\pi)^{2}}(4\pi e^{-\gamma_{E}}\mu^{2}\big{)}^{\epsilon}\,.$
(114)

Here the propagators $\rho^{\rm II}_{g}$ and $\rho^{\rm II}_{1,2}$ are
graviton and matter propagators respectively defined by

$\displaystyle\rho^{\rm II,1}_{g}$ $\displaystyle=\ell^{2}\,,$ (115)
$\displaystyle\rho^{\rm II,2}_{g}$ $\displaystyle=\big{(}\ell-
k_{23}\big{)}^{2}=\big{(}\ell+q\big{)}^{2}\,,$ $\displaystyle\rho^{\rm
II,3}_{1}$ $\displaystyle=\big{(}\ell-
k_{123}\big{)}^{2}+m_{1}^{2}=\big{(}\ell+\bar{m}_{1}u_{1}+\frac{q}{2}\big{)}^{2}+m_{1}^{2}\,,$
$\displaystyle\rho^{\rm II,4}_{2}$ $\displaystyle=\big{(}\ell-
k_{3}\big{)}^{2}+m_{2}^{2}=\big{(}\ell-\bar{m}_{2}u_{2}+\frac{q}{2}\big{)}^{2}+m_{2}^{2}\,.$

We omitted the $-i\varepsilon$ factors in the propagators.

Similarly, the integral family for the $\overline{\rm II}$ sector is given by

$F^{\overline{\rm
II}}_{i_{1},i_{2},i_{3},i_{4}}=\int_{\ell}\frac{1}{\big{(}\rho^{\overline{\rm
II},1}_{g}\big{)}^{i_{1}}\big{(}\rho^{\overline{\rm
II},2}_{g}\big{)}^{i_{2}}\big{(}\rho^{\overline{\rm
II},3}_{1}\big{)}^{i_{3}}\big{(}\rho^{\overline{\rm
II},4}_{2}\big{)}^{i_{4}}}\,,\qquad i_{1},\cdots i_{4}\in\mathbb{Z}\,,$ (116)

where the propagators are defined by

$\displaystyle\rho^{\overline{\rm II},1}_{g}$ $\displaystyle=\ell^{2}\,,$
(117) $\displaystyle\rho^{\overline{\rm II},2}_{g}$
$\displaystyle=\big{(}\ell-k_{23}\big{)}^{2}=\big{(}\ell+q\big{)}^{2}\,,$
$\displaystyle\rho^{\overline{\rm II},3}_{1}$ $\displaystyle=\big{(}\ell-
k_{123}\big{)}^{2}+m_{1}^{2}=\big{(}\ell+\bar{m}_{1}u_{1}+\frac{q}{2}\big{)}^{2}+m_{1}^{2}\,,$
$\displaystyle\rho^{\overline{\rm II},4}_{2}$ $\displaystyle=\big{(}\ell-
k_{2}\big{)}^{2}+m_{2}^{2}=\big{(}\ell+\bar{m}_{2}u_{2}+\frac{q}{2}\big{)}^{2}+m_{2}^{2}\,.$

One can check that the integral family of the $\overline{\rm II}$-sector is
related to the II-sector by crossing relation (107).

The loop integrals for the II and the $\overline{\rm II}$ diagrams are
represented by

$F^{\rm II}_{1,1,1,1}\,,\qquad F^{\overline{\rm II}}_{1,1,1,1}\,.$ (118)

The triangle diagrams can be generated from the box diagram by pinching the
box (II) diagram. The corresponding loop integrals are given by

$\displaystyle F^{\rm II}_{1,1,0,1}=F^{\overline{\rm II}}_{1,1,0,1}\,,\qquad
F^{\rm II}_{1,1,1,0}=F^{\overline{\rm II}}_{1,1,1,0}\,.$ (119)

###  4.4 Two loops

Class | couplings | power-countings
---|---|---
$\Box\\!\Box$ | $\kappa_{1}^{3}\kappa_{2}^{3}\kappa_{g}^{0}$ | $\kappa_{1}^{3}\kappa_{2}^{3}\kappa_{g}^{0}n_{1}^{2}n_{2}^{2}n_{g}^{3}+\kappa_{1}^{3}\kappa_{2}^{3}\kappa_{g}^{0}n_{1}^{2}n_{2}^{1}n_{g}^{3}+\kappa_{1}^{3}\kappa_{2}^{3}\kappa_{g}^{0}n_{1}^{1}n_{2}^{2}n_{g}^{3}$
|  | $+\kappa_{1}^{3}\kappa_{2}^{3}\kappa_{g}^{0}n_{1}^{1}n_{2}^{1}n_{g}^{3}+\kappa_{1}^{3}\kappa_{2}^{3}\kappa_{g}^{0}n_{1}^{0}n_{2}^{2}n_{g}^{3}+\kappa_{1}^{3}\kappa_{2}^{3}\kappa_{g}^{0}n_{1}^{2}n_{2}^{0}n_{g}^{3}$
$\mathbf{H}$ | $\kappa_{1}^{2}\kappa_{2}^{2}\kappa_{g}^{2}$ | $\kappa_{1}^{2}\kappa_{2}^{2}\kappa_{g}^{2}n_{1}^{1}n_{2}^{1}n_{g}^{5}+\kappa_{1}^{2}\kappa_{2}^{2}\kappa_{g}^{2}n_{1}^{1}n_{2}^{1}n_{g}^{4}$
$\mathbf{IY}$ | $\kappa_{1}^{3}\kappa_{2}^{2}\kappa_{g}^{1}$ | $\kappa_{1}^{3}\kappa_{2}^{2}\kappa_{g}^{1}n_{1}^{2}n_{2}^{1}n_{g}^{4}+\kappa_{1}^{3}\kappa_{2}^{2}\kappa_{g}^{1}n_{1}^{1}n_{2}^{1}n_{g}^{4}+\kappa_{1}^{3}\kappa_{2}^{2}\kappa_{g}^{1}n_{1}^{2}n_{2}^{0}n_{g}^{4}$
IY | $\kappa_{1}^{2}\kappa_{2}^{3}\kappa_{g}^{1}$ | $\kappa_{1}^{2}\kappa_{2}^{3}\kappa_{g}^{1}n_{1}^{1}n_{2}^{2}n_{g}^{4}+\kappa_{1}^{2}\kappa_{2}^{3}\kappa_{g}^{1}n_{1}^{1}n_{2}^{1}n_{g}^{4}+\kappa_{1}^{2}\kappa_{2}^{3}\kappa_{g}^{1}n_{1}^{0}n_{2}^{2}n_{g}^{4}$
Table 1: Power countings for classical diagrams at two loops

Let us turn to our main interest, the classification of the classical two-loop
diagrams. At the two-loop level, there exist 14 distinct power-countings,
along with four different kinds of coupling constants. According to the
relation between the classes and power countings of couplings (105), we are
led to define four classes at two loops, which are denoted as follows:
$\Box\\!\Box$, $\mathbf{H}$, $\mathbf{IY}$ and  IY . The results are
summarized in Table 1.

We define primary diagrams and sectors for each class. We start from
$\Box\\!\Box$-class. In this class, there are 6 primary diagrams, ${\rm
III}\,,{\rm IX}\,,{\rm XI}$ and their crossed diagrams $\overline{\rm
III}\,,\overline{\rm IX}\,,\overline{\rm XI}$ as depicted in Figure 8. Since
each primary diagram defines a sector, $\Box\\!\Box$-class consists of six
sectors

$\Box\\!\Box={\rm III}+{\rm IX}+{\rm XI}+\overline{\rm III}+\overline{\rm
IX}+\overline{\rm XI}\,,$ (120)

(III)

(IX)

(XI)

$(\overline{\rm III})$

$(\overline{\rm IX})$

$(\overline{\rm XI})$

Figure 8: Primary diagrams in $\Box\\!\Box$-class. These define their own
sectors.

Note that the IX and XI ($\overline{\rm IX}$ and $\overline{\rm XI}$) are a
pair of the horizontal mirror dual. Since their loop integrals are identical,
we thus omit them.

For each sector in $\Box\\!\Box$-class, an integral family is associated as
follows:

$\displaystyle F^{\rm III}_{i_{1},i_{2},\cdots i_{9}}$
$\displaystyle=\int_{\ell_{1},\ell_{2}}\frac{\big{(}\rho^{\rm
III,8}_{g}\big{)}^{-i_{8}}\big{(}\rho^{\rm
III,9}_{g}\big{)}^{-i_{9}}}{\big{(}\rho^{\rm
III,1}_{g}\big{)}^{i_{1}}\big{(}\rho^{\rm
III,2}_{g}\big{)}^{i_{2}}\big{(}\rho^{\rm
III,3}_{g}\big{)}^{i_{3}}\big{(}\rho^{\rm
III,4}_{1}\big{)}^{i_{4}}\big{(}\rho^{\rm
III,5}_{1}\big{)}^{i_{5}}\big{(}\rho^{\rm
III,6}_{2}\big{)}^{i_{6}}\big{(}\rho^{\rm III,7}_{2}\big{)}^{i_{7}}}\,,$ (121)

where

$\displaystyle\rho^{\rm III,1}_{g}$ $\displaystyle=\ell_{1}^{2}\,,$
$\displaystyle\quad\rho^{\rm III,2}_{g}$ $\displaystyle=\ell_{2}^{2}\,,$
$\displaystyle\quad\rho^{\rm III,3}_{g}$
$\displaystyle=(\ell_{1}+\ell_{2}+q)^{2}\,,$ (122) $\displaystyle\rho^{\rm
III,4}_{1}$ $\displaystyle=(\ell_{1}-k_{123})^{2}+m^{2}_{1}\,,$
$\displaystyle\quad\rho^{\rm III,5}_{1}$
$\displaystyle=(\ell_{1}+k_{1})^{2}+m^{2}_{1}\,,$ $\displaystyle\quad\rho^{\rm
III,6}_{2}$ $\displaystyle=(\ell_{1}-k_{3})^{2}+m^{2}_{2}\,,$
$\displaystyle\rho^{\rm III,7}_{2}$
$\displaystyle=(\ell_{1}-k_{2})^{2}+m^{2}_{2}\,,$
$\displaystyle\qquad\rho^{\rm III,8}_{g}$ $\displaystyle=(\ell_{1}+q)^{2},$
$\displaystyle\qquad\rho^{\rm III,9}_{g}$ $\displaystyle=(\ell_{2}+q)^{2}\;.$

and

$F^{\rm IX}_{i_{1},i_{2},\cdots
i_{9}}=\int_{\ell_{1},\ell_{2}}\frac{\big{(}\rho^{\rm
IX,8}_{g}\big{)}^{-i_{8}}\big{(}\rho^{\rm
IX,9}_{g}\big{)}^{-i_{9}}}{\big{(}\rho^{\rm
IX,1}_{g}\big{)}^{i_{1}}\big{(}\rho^{\rm
IX,2}_{g}\big{)}^{i_{2}}\big{(}\rho^{\rm
IX,3}_{g}\big{)}^{i_{3}}\big{(}\rho^{\rm
IX,4}_{1}\big{)}^{i_{4}}\big{(}\rho^{\rm
IX,5}_{1}\big{)}^{i_{5}}\big{(}\rho^{\rm
IX,6}_{2}\big{)}^{i_{6}}\big{(}\rho^{\rm IX,7}_{2}\big{)}^{i_{7}}}\,,$ (123)

where

$\displaystyle\rho^{\rm IX,1}_{g}$ $\displaystyle=\ell_{1}^{2}\,,$
$\displaystyle\quad\rho^{\rm IX,2}_{g}$ $\displaystyle=\ell_{2}^{2}\,,$
$\displaystyle\quad\rho^{\rm IX,3}_{g}$
$\displaystyle=(\ell_{1}+\ell_{2}+q)^{2},$ (124) $\displaystyle\rho^{\rm
IX,4}_{1}$ $\displaystyle=(\ell_{1}-k_{123})^{2}+m_{1}^{2}\,,$
$\displaystyle\quad\rho^{\rm IX,5}_{1}$
$\displaystyle=(\ell_{2}+k_{1})^{2}+m_{1}^{2}\,,$ $\displaystyle\quad\rho^{\rm
IX,6}_{2}$ $\displaystyle=(\ell_{2}-k_{2})^{2}+m_{2}^{2}\,,$
$\displaystyle\rho^{\rm IX,7}_{2}$
$\displaystyle=(\ell_{1}+\ell_{2}-k_{2})^{2}+m_{2}^{2}\,,$
$\displaystyle\quad\rho^{\rm IX,8}_{g}$ $\displaystyle=(\ell_{1}+q)^{2}\,,$
$\displaystyle\quad\rho^{\rm IX,9}_{g}$ $\displaystyle=(\ell_{2}+q)^{2}\,,$

Here the $\rho^{\rm III,8}_{g},\rho^{\rm III,9}_{g},\rho^{\rm IX,8}_{g}$ and
$\rho^{\rm IX,9}_{g}$ are auxiliary propagators that are not included in the
integrands of the primary diagrams. These are crucial for consistent IBP
reduction. The integral families of the crossed diagrams, $F^{\overline{\rm
III}}_{n_{1},n_{2},\cdots n_{9}}$ and $F^{\overline{\rm
IX}}_{n_{1},n_{2},\cdots n_{9}}$, can be simply derived by the crossing
relation defined in (107).

The $\mathbf{H}$-class consists of two sectors, H and $\overline{\rm H}$

$\mathbf{H}={\rm H}+\overline{\rm H}\,.$ (125)

The H and the $\overline{\rm H}$ sectors are in crossing relation and depicted
by

(H) $(\overline{\rm H})$ (126)

Their integral families $F^{\rm H}_{i}$ and $F^{\overline{\rm H}}_{i}$ are
defined by

$\displaystyle F^{\rm H}_{i_{1},\cdots,i_{9}}$
$\displaystyle=\int_{\ell_{1},\ell_{2}}\frac{\big{(}\rho^{\rm
H,8}_{1}\big{)}^{-i_{8}}\big{(}\rho^{\rm
H,9}_{2}\big{)}^{-i_{9}}}{\big{(}\rho^{\rm
H,1}_{g}\big{)}^{i_{1}}\big{(}\rho^{\rm
H,2}_{g}\big{)}^{i_{2}}\big{(}\rho^{\rm
H,3}_{g}\big{)}^{i_{3}}\big{(}\rho^{\rm
H,4}_{g}\big{)}^{i_{4}}\big{(}\rho^{\rm
H,5}_{g}\big{)}^{i_{5}}\big{(}\rho^{\rm
H,6}_{1}\big{)}^{i_{6}}\big{(}\rho^{\rm H,7}_{2}\big{)}^{i_{7}}}\,,$ (127)

where

$\displaystyle\rho^{\rm H,1}_{g}$ $\displaystyle=(\ell_{1}+q)^{2}\,,$
$\displaystyle\qquad\rho^{\rm H,2}_{g}$ $\displaystyle=(\ell_{2}+q)^{2}\,,$
$\displaystyle\qquad\rho^{\rm H,3}_{g}$ $\displaystyle=\ell_{1}^{2}\,,$ (128)
$\displaystyle\rho^{\rm H,4}_{g}$ $\displaystyle=\ell_{2}^{2}\,,$
$\displaystyle\qquad\rho^{\rm H,5}_{g}$
$\displaystyle=(\ell_{1}+\ell_{2}+q)^{2}\,,$ $\displaystyle\qquad\rho^{\rm
H,6}_{1}$ $\displaystyle=(\ell_{1}-k_{123})^{2}+m_{1}^{2}\,,$
$\displaystyle\rho^{\rm H,7}_{2}$
$\displaystyle=(\ell_{2}-k_{2})^{2}+m_{2}^{2}\,,$
$\displaystyle\qquad\rho^{\rm H,8}_{1}$
$\displaystyle=(\ell_{1}+k_{1})^{2}+m^{2}_{1}\,,$
$\displaystyle\qquad\rho^{\rm H,9}_{2}$
$\displaystyle=(\ell_{1}-k_{3})^{2}+m^{2}_{2}\,.$

The integral family of $\overline{\rm H}$ can be derived by the crossing
relation.

The $\mathbf{IY}$-class consists of 6 integral families: $\rm IY$,
$/\\!\\!\\!\rm{Y}$, $\rm YI$, $\overline{\rm IY}$,
$\overline{/\\!\\!\\!\rm{Y}}$ and $\overline{\rm YI}$, which are depicted as
follows:

(IY) (/ Y) (YI) ($\overline{\rm IY}$) ($\overline{\rm/\\!\\!\\!Y}$)
($\overline{\rm YI}$) (129)

The IY and YI ($\overline{\rm IY}$ and $\overline{\rm YI}$) sectors are
related by the horizontal mirror dual. Thus, their integral values are the
same, and we shall omit the YI and $\overline{\rm YI}$ sectors. The integral
families for the IY and the $\rm/\\!\\!\\!Y$ sectors are given by as follows:

$\displaystyle F^{\rm IY}_{i_{1},i_{2},\cdots,i_{9}}$
$\displaystyle=\int_{\ell_{1},\ell_{2}}\frac{\big{(}\rho^{\rm
IY,8}_{2}\big{)}^{-i_{8}}\big{(}\rho^{\rm
IY,9}_{g}\big{)}^{-i_{9}}}{\big{(}\rho^{\rm
IY,1}_{g}\big{)}^{i_{1}}\big{(}\rho^{\rm
IY,2}_{g}\big{)}^{i_{2}}\big{(}\rho^{\rm
IY,3}_{g}\big{)}^{i_{3}}\big{(}\rho^{\rm
IY,4}_{g}\big{)}^{i_{4}}\big{(}\rho^{\rm
IY,5}_{1}\big{)}^{i_{5}}\big{(}\rho^{\rm
IY,6}_{1}\big{)}^{i_{6}}\big{(}\rho^{IY,7}_{2}\big{)}^{i_{7}}}\,,$ (130)

where

$\displaystyle\rho^{\rm IY,1}_{g}$ $\displaystyle=(\ell_{1}+q)^{2}\,,$
$\displaystyle\quad\rho^{\rm IY,2}_{g}$ $\displaystyle=\ell_{1}^{2}\,,$
$\displaystyle\quad\rho^{\rm IY,3}_{g}$ $\displaystyle=\ell_{2}^{2}\,,$ (131)
$\displaystyle\rho^{\rm IY,4}_{g}$
$\displaystyle=(\ell_{1}+\ell_{2}+q)^{2}\,,$
$\displaystyle\quad\tilde{\rho}^{\rm IY,5}_{1}$
$\displaystyle=(\ell_{2}-k_{123})^{2}+m_{1}^{2}\,,$
$\displaystyle\quad\tilde{\rho}^{\rm IY,6}_{1}$
$\displaystyle=(\ell_{1}+k_{1})^{2}+m_{1}^{2}\,,$
$\displaystyle\tilde{\rho}^{\rm IY,7}_{2}$
$\displaystyle=(\ell_{1}-k_{2})^{2}+m_{2}^{2}\,,$
$\displaystyle\quad\tilde{\rho}^{\rm IY,8}_{2}$
$\displaystyle=(\ell_{2}-k_{3})^{2}+m_{2}^{2}\;,$ $\displaystyle\quad\rho^{\rm
IY,9}_{g}$ $\displaystyle=(\ell_{2}+q)^{2}\,.$

and

$F^{\rm/\\!\\!\\!Y}_{i_{1},i_{2},\cdots,i_{9}}=\int_{\ell_{1},\ell_{2}}\frac{\big{(}\rho^{\rm/\\!\\!\\!Y,8}_{2}\big{)}^{-i_{8}}\big{(}\rho^{\rm/\\!\\!\\!Y,9}_{g}\big{)}^{-i_{9}}}{\big{(}\rho^{\rm/\\!\\!\\!Y,1}_{g}\big{)}^{i_{1}}\big{(}\rho^{\rm/\\!\\!\\!Y,2}_{g}\big{)}^{i_{2}}\big{(}\rho^{\rm/\\!\\!\\!Y,3}_{g}\big{)}^{i_{3}}\big{(}\rho^{\rm/\\!\\!\\!Y,4}_{g}\big{)}^{i_{4}}\big{(}\rho^{\rm/\\!\\!\\!Y,5}_{1}\big{)}^{i_{5}}\big{(}\tilde{\rho}^{\rm/\\!\\!\\!Y,6}_{1}\big{)}^{i_{6}}\big{(}\rho^{\rm/\\!\\!\\!Y,7}_{2}\big{)}^{i_{7}}}\,,$
(132)

where

$\displaystyle\rho^{\text{/\\!\\!\\!\rm{Y}},1}_{g}$
$\displaystyle=(\ell_{1}+q)^{2}\,,$
$\displaystyle\quad\rho^{\text{/\\!\\!\\!\rm{Y}},2}_{g}$
$\displaystyle=\ell_{1}^{2}\,,$
$\displaystyle\quad\rho^{\text{/\\!\\!\\!\rm{Y}},3}_{g}$
$\displaystyle=\ell_{2}^{2}\,,$ (133)
$\displaystyle\rho^{\text{/\\!\\!\\!\rm{Y}},4}_{g}$
$\displaystyle=(\ell_{1}+\ell_{2}+q)^{2}\,,$
$\displaystyle\quad\rho^{\text{/\\!\\!\\!\rm{Y}},5}_{1}$
$\displaystyle=(\ell_{2}-k_{123})^{2}+m_{1}^{2}\,,$
$\displaystyle\quad\rho^{\text{/\\!\\!\\!\rm{Y}},6}_{1}$
$\displaystyle=(\ell_{1}+\ell_{2}-k_{123})^{2}+m_{1}^{2}\,,$
$\displaystyle\rho^{\text{/\\!\\!\\!\rm{Y}},7}_{2}$
$\displaystyle=(\ell_{1}-k_{2})^{2}+m_{2}^{2}\,,$
$\displaystyle\quad\rho^{\text{/\\!\\!\\!\rm{Y}},8}_{2}$
$\displaystyle=(\ell_{1}-k_{3})^{2}+m_{2}^{2}\;,$
$\displaystyle\quad\rho^{\text{/\\!\\!\\!\rm{Y}},9}_{g}$
$\displaystyle=(\ell_{2}+q)^{2}\,.$

Finally, the  IY -class consists of six sectors,
$\backslash\\!\\!\\!\\!\\!\scalebox{-1.0}[1.0]{\rotatebox[origin={c}]{180.0}{
\rm{Y}}}$,  YI ,
$\overline{\scalebox{-1.0}[1.0]{\rotatebox[origin={c}]{180.0}{\rm IY}}}$,
$\overline{\backslash\\!\\!\\!\scalebox{-1.0}[1.0]{\rotatebox[origin={c}]{180.0}{\rm{Y}}}}$
and $\overline{\scalebox{-1.0}[1.0]{\rotatebox[origin={c}]{180.0}{\rm YI}}}$,
which are illustrated as follows:

(  IY ) (  / Y )
$(\scalebox{-1.0}[1.0]{\rotatebox[origin={c}]{180.0}{YI}})$
($\overline{\scalebox{-1.0}[1.0]{\rotatebox[origin={c}]{180.0}{\rm IY}}}$)
($\overline{\backslash\\!\\!\\!\scalebox{-1.0}[1.0]{\rotatebox[origin={c}]{180.0}{Y}}}$)
($\overline{\scalebox{-1.0}[1.0]{\rotatebox[origin={c}]{180.0}{\rm YI}}}$)
(134)

One can easily notice that the  IY -class is the vertical mirror image of the
$\mathbf{IY}$-class. This implies that the loop integrals of  IY -class can be
derived by the mapping for the vertical mirror dual in (108) without explicit
computation from the results of $\mathbf{IY}$-class.

##  5 Soft Expansion and IBP reduction

In the previous section, we introduced the new framework for separating the
classical diagrams and their classification. Nonetheless, one must also
account for the residual quantum contributions present in the numerators of
loop integrands, even for the classical diagrams. These quantum contributions
arise from the terms with higher powers of momentum transfer $q^{\mu}$ or loop
momenta $\ell_{i}$ and do not contribute to the classical long-range
interaction – these correspond to UV properties of GR Bern:2019crd ;
Bern:2019nnu .

This section briefly reviews the soft expansion and IBP reduction framework
established in Parra-Martinez:2020dzs  in the context of the off-shell
recursions. We also apply the framework to the one-loop classical amplitude.

###  5.1 Brief review on the soft expansion

To obtain the classical contributions for a given loop integrand for a
classical diagram, one typically employs the soft expansion in the context of
the method of regions Beneke:1997zp . There are four different regions for the
internal graviton loop momenta: hard, soft, potential and radiation. The
potential region corresponds to the conservative dynamics. In this region the
loop momenta $\ell=(\omega,\vec{\ell})$ scales as

$(\omega,\vec{\ell})\sim(|\vec{q}||\vec{v}|,|\vec{q}|)\,.$ (135)

The scales in the soft region satisfy the hierarchy
$m_{1},m_{2},\bar{k}_{i}\gg q^{\mu},\ell_{i}$. We implement a scaling
transformation on the variables in $\hbar$ as

$\displaystyle\ell_{i}\to\hbar\ell_{i}\,,\qquad q\to\hbar q\,,\qquad u_{i}\to
u_{i}\,,\qquad m_{i}\to m_{i}\,.$ (136)

Upon this scaling, a loop integral $I$ for a classical diagram is expanded in
$\hbar$

$I=\cdots+\hbar^{-2}I^{(-2)}+\hbar^{-1}I^{(-1)}+\hbar^{0}I^{(0)}+\sum_{n>0}\hbar^{n}I^{(n)}\,.$
(137)

Obviously, the terms with positive powers in $\hbar$ are the quantum
contributions, and we should truncate them. The remaining $\hbar^{0}$ and
$\hbar^{n<0}$ order terms are the so-called classical and superclassical
terms, respectively.

Let us consider in more detail the scaling of the integrands. Each integration
measure scales as $\hbar^{4}$ in $4$-dimensions

$\int_{\ell_{i}}\to\hbar^{4}\int_{\ell_{i}}\,.$ (138)

Regarding the propagators, the graviton propagators scale homogeneously as
$\hbar^{-2}$

$\frac{1}{\ell_{i}^{2}}\to\frac{\hbar^{-2}}{\ell_{i}^{2}}\,,\qquad\frac{1}{(\ell_{i}\pm
q)^{2}}\to\frac{\hbar^{-2}}{(\ell_{i}\pm q)^{2}}\,.$ (139)

On the other hand, the matter propagators scale inhomogeneously and can be
expanded in $\hbar$ as

$\displaystyle\frac{1}{\ell^{2}+2\ell\cdot(\frac{q}{2}\pm\bar{m}_{I}u_{I})}\to\frac{\pm
1}{2\bar{m}_{I}(u_{I}\cdot\ell)}\left[\hbar^{-1}\mp\frac{\ell^{2}+\ell\cdot
q}{2\bar{m}_{I}(\ell\cdot u_{I})}+\hbar\left(\frac{\ell^{2}+\ell\cdot
q}{2\bar{m}_{I}(\ell\cdot u_{I})}\right)^{2}+\cdots\right]\,,$ (140)

and the scaling of the matter propagators start from the $\hbar^{-1}$-order.

The expansion of the matter propagator is conveniently represented by
introducing a linear matter propagator, which scales precisely as
$\hbar^{-1}$,

$\frac{\pm 1}{2\bar{m}_{I}(u_{I}\cdot\ell)}\,,$ (141)

and the remainder in the expansion (140) is absorbed into the numerator. Let
us consider a loop integrand $I$ for a classical diagram with a numerator
$\mathcal{N}$

$I=\int_{\ell_{1},\cdots,\ell_{l}}\frac{\mathcal{N}[\ell_{1},\cdots,\ell_{l};u_{1},u_{2},q^{\mu}]}{\big{(}\rho^{1}_{g}\big{)}^{i_{1}}\cdots\big{(}\rho^{m}_{g}\big{)}^{i_{m}}\big{(}\rho^{m+1}_{1}\big{)}^{i_{m+1}}\cdots\big{(}\rho_{1}^{n}\big{)}^{i_{n}}\big{(}\rho_{2}^{i_{n+1}}\big{)}^{i_{n+1}}\cdots\big{(}\rho_{2}^{p}\big{)}^{i_{p}}}\,,$
(142)

where $1<m<n<p$ and $\rho^{i}_{g}$, and $\rho^{i}_{1,2}$ are the gravitons and
the full matter propagators respectively.

After the soft expansion on the loop integrands, we may recast the integrand
using the modified numerator $\tilde{\mathcal{N}}$, which includes all the
terms arising from the expansion of the matter propagators (140)

$I=\int_{\ell_{1},\cdots,\ell_{l}}\frac{\tilde{\mathcal{N}}[\ell_{1},\cdots,\ell_{l};u_{1},u_{2},q^{\mu}]}{\big{(}\rho^{1}_{g}\big{)}^{i_{1}}\cdots\big{(}\rho^{m}_{g}\big{)}^{i_{m}}\big{(}\tilde{\rho}^{m+1}_{1}\big{)}^{i_{m+1}}\cdots\big{(}\tilde{\rho}_{1}^{n}\big{)}^{i_{n}}\big{(}\tilde{\rho}_{2}^{i_{n+1}}\big{)}^{i_{n+1}}\cdots\big{(}\tilde{\rho}_{2}^{p}\big{)}^{i_{p}}}\,,$
(143)

where $\tilde{\rho}^{i}_{1}$ and $\tilde{\rho}^{i}_{2}$ are the linearized
matter propagators defined in (141). To get a consistent classical
contribution, we truncate the numerators $\tilde{\mathcal{N}}$ by imposing the
power counting constraint after the $\hbar$ expansion (136). If
$\tilde{\mathcal{N}}$ scales as
$\tilde{\mathcal{N}}\to\sum_{N}\hbar^{N}\tilde{\mathcal{N}}^{{\scriptscriptstyle(N)}}$,
we keep the terms satisfying the following relation to truncate the quantum
part:

$4l-2n_{g}-n_{1}-n_{2}+N\leq 0\,,$ (144)

where $n_{g}$ and $n_{1,2}$ are the number of graviton and the scalar fields
respectively. Accordingly, we have $\tilde{\mathcal{N}}_{\rm classical}$ after
the truncation as

$\tilde{\mathcal{N}}_{\rm classical}=\sum_{N\leq
2n_{g}+n_{1}+n_{2}-4l}\hbar^{N}\tilde{\mathcal{N}}^{{\scriptscriptstyle(N)}}\,.$
(145)

As the matter propagators are replaced with the linear ones, we have to modify
the integral families in (103) by accommodating the linearized matter
propagators, which we shall call the soft integral family,

$\displaystyle
j_{i_{1},\cdots,i_{p}}=\int_{\ell_{1},\ell_{2},\cdots,\ell_{l}}\frac{1}{\big{(}\rho^{1}_{g}\big{)}^{i_{1}}\cdots\big{(}\rho^{m}_{g}\big{)}^{i_{m}}\big{(}\tilde{\rho}^{m+1}_{1}\big{)}^{i_{m+1}}\cdots\big{(}\tilde{\rho}_{1}^{n}\big{)}^{i_{n}}\big{(}\tilde{\rho}_{2}^{i_{n+1}}\big{)}^{i_{n+1}}\cdots\big{(}\tilde{\rho}_{2}^{p}\big{)}^{i_{p}}}\,.$
(146)

Nonetheless, the intrinsic relationship between the sectors and the integral
families remains intact even after the soft expansion. Accordingly, we will
exploit IBP reduction with respect to the soft integral families.

The IBP reduction allows us to represent a loop integrand as a linear
combination of a set of basis integrals called master integrals. A complete
set of master integrals may generate a system of differential equations by
acting on the following differential operator to the master integrals:

$\frac{\mathrm{d}}{\mathrm{d}y}=\frac{2x^{2}}{x^{2}-1}\frac{\mathrm{d}}{\mathrm{d}x}\,.$
(147)

Here $x$ and $y$ are defined in (93) and (95), and master integrals are
functions of these variables. We can evaluate the loop integrals by solving
the differential equations with a proper boundary condition. We may adopt the
results in Parra-Martinez:2020dzs .

###  5.2 One loop example

Let us consider the one-loop example for elucidating our approach. As
discussed in Section 4.3.3, the only contribution to the classical amplitude
at one loop is from the box and triangle diagrams. We shall collect the
integrands associated with the box and triangle diagrams using the power
counting rules in (109) and (110). At one loop level, there is only $\Box$
class, which is decomposed into two sectors, the II and $\overline{\rm II}$
sectors defined in (111). The soft integral families for each sector, (112)
and (116), are given by

$\displaystyle j^{\rm II}_{i_{1},i_{2},i_{3},i_{4}}$
$\displaystyle=\int_{\ell}\frac{1}{\mathcal{D}^{\rm
II}_{i_{1},i_{2},i_{3},i_{4}}}\,,\qquad\mathcal{D}^{\rm
II}_{i_{1},i_{2},i_{3},i_{4}}=\big{(}{\rho}^{\rm
II,1}_{g}\big{)}^{i_{1}}\big{(}{\rho}^{\rm
II,2}_{g}\big{)}^{i_{2}}\big{(}\tilde{\rho}^{\rm
II,3}_{1}\big{)}^{i_{3}}\big{(}\tilde{\rho}^{\rm II,4}_{2}\big{)}^{i_{4}}\,,$
(148) $\displaystyle j^{\overline{\rm II}}_{i_{1},i_{2},i_{3},i_{4}}$
$\displaystyle=\int_{\ell}\frac{1}{\mathcal{D}^{\overline{\rm
II}}_{i_{1},i_{2},i_{3},i_{4}}}\,,\qquad\mathcal{D}^{\overline{\rm
II}}_{i_{1},i_{2},i_{3},i_{4}}=\big{(}{\rho}^{\overline{\rm
II},1}_{g}\big{)}^{i_{1}}\big{(}{\rho}^{\overline{\rm
II},2}_{g}\big{)}^{i_{2}}\big{(}\tilde{\rho}^{\overline{\rm
II},3}_{1}\big{)}^{i_{3}}\big{(}\tilde{\rho}^{\overline{\rm
II},4}_{2}\big{)}^{i_{4}}\,,$

where $\rho^{{\rm II}}_{g}$ and ${\rho}^{\overline{\rm II}}_{g}$ are the
graviton propagators in each sector, and $\tilde{\rho}^{{\rm II}}_{1}$ and
$\tilde{\rho}^{{\rm II}}_{2}$ ($\tilde{\rho}^{\overline{\rm II}}_{1}$ and
$\tilde{\rho}^{\overline{\rm II}}_{2}$) are the linear matter propagators in
the II sector ($\overline{\rm II}$ sector) for the scalar field 1 and 2
respectively

$\displaystyle\tilde{\rho}^{\rm II,3}_{1}$ $\displaystyle=2\ell\cdot u_{1}\,,$
$\displaystyle\qquad\tilde{\rho}^{\rm II,4}_{2}$ $\displaystyle=-2\ell\cdot
u_{2}\,,$ (149) $\displaystyle\tilde{\rho}^{\overline{\rm II},3}_{1}$
$\displaystyle=2\ell\cdot u_{1}\,,$
$\displaystyle\qquad\tilde{\rho}^{\overline{\rm II},4}_{2}$
$\displaystyle=2\ell\cdot u_{2}\,.$

One can check that the integral families (148) are related to each other by
the crossing relation (107).

####  5.2.1 II sector

We first consider the II sector. From the IBP reduction of the soft integral
family, there are three master integrals $f^{\rm II}_{1}$, $f^{\rm II}_{2}$
and $f^{\rm II}_{3}$ defined as

$\displaystyle f^{\rm II}_{1}=\epsilon q^{2}j_{0,0,2,1}\,,\qquad f^{\rm
II}_{2}=-\epsilon|q|j_{1,0,1,1}\,,\qquad f^{\rm
II}_{3}=\epsilon^{2}\sqrt{y^{2}-1}q^{2}j_{1,1,1,1}\,,$ (150)

where $\epsilon$ is the dimensional regularization parameter, $d=4-2\epsilon$.
Acting on the differential operator (147), we derive a differential equation

$\frac{\mathrm{d}\vec{f}_{\rm
II}}{\mathrm{~{}d}x}=\epsilon\frac{A}{x}\vec{f}_{\rm II}\,,\qquad\vec{f}_{\rm
II}=\big{\\{}f^{\rm II}_{1}\,,f^{\rm II}_{2}\,,f^{\rm II}_{3}\big{\\}}\,,$
(151)

where $A$ is a $3\times 3$ constant matrix

$A=\begin{pmatrix}0&0&0\\\ 0&0&0\\\ 1&0&0\end{pmatrix}\,.$ (152)

In Parra-Martinez:2020dzs ; Bern:2019crd ; Beenakker:1988jr , the solution of
the differential equation (151) is constructed in the potential region with a
boundary condition at the static limit $y\to 1$

$\displaystyle f_{1}^{\rm II}$ $\displaystyle=0\,,$ (153) $\displaystyle
f_{2}^{\rm II}$
$\displaystyle=\epsilon^{2}\left(q^{2}\right)^{-\epsilon}e^{\gamma_{\mathrm{E}}\epsilon}\frac{\sqrt{\pi}\Gamma\left(\frac{1}{2}-\epsilon\right)^{2}\Gamma\left(\epsilon+\frac{1}{2}\right)}{2\Gamma(1-2\epsilon)}\,,$
$\displaystyle f_{3}^{\rm II}$
$\displaystyle=\epsilon^{2}\left(q^{2}\right)^{-\epsilon}e^{\gamma_{\mathrm{E}}\epsilon}\frac{\mathrm{i}\pi}{2}\frac{\Gamma(-\epsilon)^{2}\Gamma(1+\epsilon)}{\Gamma(-2\epsilon)}\,.$

Since $f^{\rm II}_{1}$ does not contribute to the result at the potential
region, we ignore it from now on.

From the solution of the one-loop recursions, we can derive the explicit form
of the one-loop soft integral for the II sector, which is decomposed into
three parts:

$\displaystyle\mathcal{M}^{\rm II}$
$\displaystyle=\int_{\ell}\frac{\tilde{\mathcal{N}}^{\rm
II}_{1,1,1,1}}{\mathcal{D}^{\rm
II}_{1,1,1,1}}+\int_{\ell}\frac{\tilde{\mathcal{N}}^{\rm
II}_{1,1,0,1}}{\mathcal{D}^{\rm
II}_{1,1,0,1}}+\int_{\ell}\frac{\tilde{\mathcal{N}}^{\rm
II}_{1,1,1,0}}{\mathcal{D}^{\rm II}_{1,1,1,0}}\,,$ (154)

where $\mathcal{D}^{\rm II}_{i_{1},i_{2},i_{3},i_{4}}$ is the propagators
defined in (148). The first term corresponds to the II diagram, and the second
and third correspond to the $\triangle$ and  $\triangle$ diagrams. The
numerators $\tilde{\mathcal{N}}^{\rm II}$ are given by

$\displaystyle\tilde{\mathcal{N}}^{\rm II}_{1,1,1,1}$
$\displaystyle=\frac{4m_{1}^{4}m_{2}^{4}\left(2\sigma^{2}-1\right)^{2}}{\bar{m}_{1}\bar{m}_{2}}$
(155)
$\displaystyle\quad+2m_{1}^{2}m_{2}^{2}\left(2\sigma^{2}-1\right)\Big{(}\frac{m_{1}}{\ell\cdot
u_{2}}-\frac{m_{2}}{\ell\cdot u_{1}}\Big{)}\left[L+8\sigma(\ell\cdot
u_{1})(\ell\cdot u_{2})\right]$
$\displaystyle\quad+m_{1}m_{2}\bigg{[}\frac{\big{(}m_{1}L-4m_{2}\sigma(\ell\cdot
u_{2})^{2}\big{)}^{2}}{\left(\ell\cdot
u_{2}\right){}^{2}}+\frac{m_{2}L\left(8m_{2}\sigma\left(\ell\cdot
u_{2}\right){}^{2}-m_{1}L\right)}{(\ell\cdot u_{1})(\ell\cdot u_{2})}$
$\displaystyle\qquad\qquad\quad+\frac{m_{2}^{2}L^{2}}{\left(\ell\cdot
u_{1}\right)^{2}}+\frac{8m_{1}\sigma(\ell\cdot
u_{1})\left(m_{1}L-4m_{2}\sigma(\ell\cdot u_{2})^{2}\right)}{\ell\cdot u_{2}}$
$\displaystyle\qquad\qquad\quad+16m_{1}^{2}\sigma^{2}\left(\ell\cdot
u_{1}\right){}^{2}+4m_{1}m_{2}q^{2}\sigma\left(1-2\sigma^{2}\right)\bigg{]}\,,$
$\displaystyle\tilde{\mathcal{N}}^{\rm II}_{1,1,0,1}$
$\displaystyle=\frac{2m_{2}^{2}m_{1}^{4}}{q^{2}\bar{m}_{1}}\bigg{[}\sigma^{2}(2\ell+q)^{2}+q^{2}\left(5\sigma^{2}-1\right)\bigg{]}$
$\displaystyle\quad-\frac{m_{1}^{2}m_{2}^{2}(\ell\cdot
u_{1})}{q^{2}}\bigg{[}16\left(\ell\cdot
u_{2}\right){}^{2}+(2\ell+q)^{2}+q^{2}\left(5-4\sigma^{2}\right)\bigg{]}$
$\displaystyle\quad-\frac{m_{1}^{2}m_{2}^{2}}{(\ell\cdot
u_{1})q^{2}}\Big{(}\ell\cdot(\ell+q)\Big{)}\bigg{[}\sigma^{2}(2\ell+q)^{2}+q^{2}\left(5\sigma^{2}-1\right)+2\left(\ell\cdot
u_{2}\right)^{2}\bigg{]}$
$\displaystyle\quad+\frac{4m_{1}^{2}m_{2}^{2}}{q^{2}}\bigg{[}4m_{1}\sigma^{2}\left(\ell\cdot
u_{1}\right)^{2}+m_{1}\left(\ell\cdot
u_{2}\right)^{2}-\sigma\Big{(}2\ell\cdot(\ell+q)+5q^{2}\Big{)}(\ell\cdot
u_{2})\bigg{]}\,,$ $\displaystyle\tilde{\mathcal{N}}^{\rm II}_{1,1,1,0}$
$\displaystyle=\mathcal{N}^{\rm II}_{1,1,0,1}\Big{|}_{u_{1}\leftrightarrow
u_{2},m_{1}\leftrightarrow m_{2},q\to-q}\,.$

where $L=\left(2\sigma^{2}-1\right)(\ell^{2}+\ell\cdot q)$.

Applying IBP reduction for tensor reduction, we can represent
$\mathcal{M}^{\rm II}$ in terms of the master integrals

$\mathcal{M}^{\rm II}=\sum_{i=2}^{3}c^{\rm II}_{i}f^{\rm II}_{i}\big{|}_{\rm
non-analytic}\,.$ (156)

We should keep only non-analytic terms in $|q|$ or $\epsilon$ since the
analytic terms are irrelevant for the classical long-range interaction Parra-
Martinez:2020dzs . The coefficients $c^{\rm II}_{i}$, which are independent on
the loop momentum $\ell$, are

$\displaystyle c^{\rm II}_{2}$
$\displaystyle=\frac{4m_{1}^{2}m_{2}^{2}\left(m_{1}{+}m_{2}\right)}{|q|}\Bigg{[}\frac{\left(2\sigma^{2}-1\right)^{2}}{\sigma-1}-\frac{2\left(2\sigma^{2}-1\right)}{\epsilon}+\frac{(15\sigma^{2}+3)+4\left(4\sigma^{2}-1\right)\epsilon}{8\epsilon(\epsilon-1)}\Bigg{]}\,,$
(157) $\displaystyle c^{\rm II}_{3}$
$\displaystyle=\frac{4m_{1}^{3}m_{2}^{3}\left(1-2\sigma^{2}\right)^{2}}{\sqrt{\sigma^{2}-1}\epsilon^{2}q^{2}}-\frac{m_{1}m_{2}\left(2\sigma^{2}-1\right)^{2}\left(2m_{2}m_{1}\sigma+m_{1}^{2}+m_{2}^{2}\right)}{\left(\sigma^{2}-1\right)^{3/2}\epsilon}\,.$

Substituting the coefficients and the master integrals into (156), we have

$\displaystyle\mathcal{M}^{\rm II}$
$\displaystyle=\frac{4\pi^{2}m_{1}^{2}m_{2}^{2}\sigma\left(2\sigma^{2}-1\right)(m_{1}+m_{2})}{(q^{2})^{\epsilon+\frac{1}{2}}}-\frac{4i\pi
m_{1}^{3}m_{2}^{3}(2\sigma^{2}-1)^{2}}{\epsilon\sqrt{\sigma^{2}-1}(q^{2})^{\epsilon+1}}$
(158)
$\displaystyle\quad+\frac{3\pi^{2}m_{1}^{2}m_{2}^{2}(m_{1}+m_{2})\left(5\sigma^{2}-1\right)}{4(q^{2})^{\frac{1}{2}+\epsilon}}\,,$

where the first line is from the box and the crossed box diagrams, and the
second line is the contribution from the triangle diagrams.

####  5.2.2 $\overline{\rm II}$ sector

We next consider the $\overline{\rm II}$ sector. According to the crossing
relation (107), the definition of the master integrals is replaced by

$f^{\overline{\rm II}}_{1}=\epsilon q^{2}j^{\overline{\rm
II}}_{0,0,2,1}\,,\qquad f^{\overline{\rm II}}_{2}=-\epsilon|q|j^{\overline{\rm
II}}_{1,0,1,1}\,,\qquad f^{\overline{\rm
II}}_{3}=-\epsilon^{2}\sqrt{y^{2}-1}q^{2}j^{\overline{\rm II}}_{1,1,1,1}\,.$
(159)

On the other hand, the differential equation is not altered under the crossing
rule. The difference arises from the boundary condition at the static limit,
$y\to 1$. Solving the differential equation, we have the master integrals
Parra-Martinez:2020dzs

$\displaystyle f_{1}^{\overline{\rm II}}$ $\displaystyle=0\,,$ (160)
$\displaystyle f_{2}^{\overline{\rm II}}$
$\displaystyle=\epsilon^{2}\left(q^{2}\right)^{-\epsilon}e^{\gamma_{\mathrm{E}}\epsilon}\frac{\sqrt{\pi}\Gamma\left(\frac{1}{2}-\epsilon\right)^{2}\Gamma\left(\epsilon+\frac{1}{2}\right)}{2\Gamma(1-2\epsilon)}\,,$
$\displaystyle f_{3}^{\overline{\rm II}}$ $\displaystyle=0\,.$

The explicit form of the one-loop integrand for the $\overline{\rm II}$ sector
is also derived by solving the one-loop recursions and divided into three
parts:

$\displaystyle\mathcal{M}^{\overline{\rm II}}$
$\displaystyle=\int_{\ell}\frac{\mathcal{N}^{\overline{\rm
II}}_{1,1,1,1}}{\mathcal{D}^{\overline{\rm
II}}_{1,1,1,1}}+\int_{\ell}\frac{\mathcal{N}^{\overline{\rm
II}}_{1,1,0,1}}{\mathcal{D}^{\overline{\rm
II}}_{1,1,0,1}}+\int_{\ell}\frac{\mathcal{N}^{\overline{\rm
II}}_{1,1,1,0}}{\mathcal{D}^{\overline{\rm II}}_{1,1,1,0}}\,,$ (161)

where $\mathcal{D}^{\overline{\rm II}}_{i_{1},i_{2},i_{3},i_{4}}$ is defined
in (148) and the numerators are given by

$\displaystyle\mathcal{N}^{\overline{\rm II}}_{1,1,1,1}$
$\displaystyle=\mathcal{N}^{{\rm II}}_{1,1,1,1}\Big{|}_{\rm crossing}\,,$
(162) $\displaystyle\mathcal{N}^{\overline{\rm II}}_{1,1,0,1}$
$\displaystyle=\mathcal{N}^{{\rm II}}_{1,1,0,1}\,,$
$\displaystyle\mathcal{N}^{\overline{\rm II}}_{1,1,1,0}$
$\displaystyle=\mathcal{N}^{{\rm II}}_{1,1,1,0}\,,$

Again, applying IBP reduction to (161) for tensor reduction,
$\mathcal{M}^{\overline{\rm II}}$ reduces to

$\mathcal{M}^{\overline{\rm II}}=c^{\overline{\rm II}}_{2}f^{\overline{\rm
II}}_{2}\Big{|}_{\rm non-analytc}\,,$ (163)

where the coefficients are

$\displaystyle c^{\overline{\rm II}}_{2}$
$\displaystyle=\frac{4m_{1}^{2}m_{2}^{2}\left(m_{1}{+}m_{2}\right)}{|q|}\Bigg{[}\frac{\left(2\sigma^{2}-1\right)^{2}}{\sigma-1}+\frac{2\left(2\sigma^{2}-1\right)}{\epsilon}+\frac{15\sigma^{2}+3}{8\epsilon(\epsilon-1)}+\frac{4\sigma^{2}-1}{2(\epsilon-1)}\Bigg{]}\,.$
(164)

By substituting the result into (163), we have the amplitude from the
$\overline{\rm II}$ sector

$\displaystyle\mathcal{M}^{\overline{\rm
II}}=-\frac{4\pi^{2}m_{1}^{2}m_{2}^{2}\sigma\left(2\sigma^{2}-1\right)(m_{1}+m_{2})}{(q^{2})^{\epsilon+\frac{1}{2}}}+\frac{3\pi^{2}m_{1}^{2}m_{2}^{2}(m_{1}+m_{2})\left(5\sigma^{2}-1\right)}{4(q^{2})^{\epsilon+\frac{1}{2}}}\,.$
(165)

####  5.2.3 Total classical one-loop amplitude

The total classical one-loop amplitude is given by combining the results from
the II and $\overline{\rm II}$ sectors,

$\mathcal{M}^{(1)}_{\rm pot}=4G^{2}\Big{(}\mathcal{M}^{\rm
II}+\mathcal{M}^{\overline{\rm II}}\Big{)}\,.$ (166)

Here the overall factor arises from the $(-i)$ factor in the LSZ reduction
formula (56), 1-loop coupling $\kappa^{2}=8\pi G$ and the normalization of the
loop integral (114), $-i\frac{i(8\pi G)^{2}}{(4\pi)^{2}}=4G^{2}$.

Adding the contributions from the II and $\overline{\rm II}$ sectors, (158)
and (165), and taking $\epsilon\to 0$ limit, we obtain total classical one-
loop amplitude in the potential region using (166)

$\displaystyle\mathcal{M}^{{\scriptscriptstyle(1)}}_{\rm pot}$
$\displaystyle=-\frac{16iG^{2}\pi
m_{1}^{3}m_{2}^{3}\left(2\sigma^{2}-1\right)^{2}}{\sqrt{\sigma^{2}-1}q^{2}}+\frac{6G^{2}\pi^{2}m_{1}^{2}m_{2}^{2}(m_{1}+m_{2})\left(5\sigma^{2}-1\right)}{|q|}\,.$
(167)

The first term on the RHS is from the box and crossed box diagrams, and the
last is from the triangles. The box diagram contribution is superclssical and
does not contribute to the conservative potential – it vanishes by the
iteration. This result is consistent with the known result Bjerrum-
Bohr:2021vuf .

##  6 Two Loop Amplitude

In this section, we calculate the classical loop amplitude at two loops,
combining all the ingredients we developed so far. According to the power-
countings of the couplings, we identified the four classes, $\Box\\!\Box$,
$\mathbf{H}$, $\mathbf{IY}$ and  IY , and the corresponding sectors for each
class in Section 4.4. We represent each loop integrand as a linear combination
of the master integrals using IBP reduction. We employed the LiteRed
Lee:2012cn ; Lee:2013mka ; Smirnov:2019qkx  based on the Laporta program
Laporta:2000dsw .

For each sector, we introduce the soft integral family and the master
integrals. We also present the differential equations for the soft integral
families and solutions at the potential region. The solutions are obtained in
Smirnov:2001cm ; Henn:2013woa ; Bianchi:2016yiq ; Parra-Martinez:2020dzs , and
we convert the result in our convention. Finally, we derive the total
classical amplitude by collecting all the amplitude contributions for each
sector as follows:

$\mathcal{M}^{{\scriptscriptstyle(2)}}_{\rm
pot}=\frac{2G^{3}}{\pi}\Big{(}\mathcal{M}^{\Box\\!\Box}+\mathcal{M}^{\mathbf{H}}+\mathcal{M}^{\mathbf{IY}}+\mathcal{M}^{\mathbf{IY}}\big{|}_{m_{1}\leftrightarrow
m_{2}}\Big{)}\,.$ (168)

###  6.1 Double box ($\Box\\!\Box$) class

Six sectors are identified within the $\Box\\!\Box$ class: III, IX, XI, and
their crossings, as described in (120). Since the IX and XI sectors are in the
horizontal mirror dual, as well as their crossed counterparts $\overline{\rm
IX}$ and $\overline{\rm XI}$, we will consider the IX and $\overline{\rm IX}$
sectors only. The overall contribution to the amplitude from the $\Box\\!\Box$
class is

$\mathcal{M}^{\Box\\!\Box}=\mathcal{M}^{\rm III}+2\mathcal{M}^{\rm
IX}+\mathcal{M}^{\overline{\rm III}}+2\mathcal{M}^{\overline{\rm IX}}\,.$
(169)

####  6.1.1 $\rm III$ and $\overline{\rm III}$ sectors

The integral family associated with the III sector is introduced in (121). The
corresponding soft integral family is obtained by replacing the full matter
propagators $\rho^{\rm III}_{1}$ and $\rho^{\rm III}_{1}$ with the linearized
propagators, $\tilde{\rho}^{\rm III}_{1}$ and $\tilde{\rho}^{\rm III}_{2}$,

$\displaystyle
j^{\text{III}}_{i_{1},i_{2},\cdots,i_{9}}=\int_{\ell_{1},\ell_{2}}\frac{\big{(}\rho^{\rm
III,8}_{g}\big{)}^{-i_{8}}\big{(}\rho^{\rm
III,9}_{g}\big{)}^{-i_{9}}}{\big{(}\rho^{\rm
III,1}_{g}\big{)}^{i_{1}}\big{(}\rho^{\rm
III,2}_{g}\big{)}^{i_{2}}\big{(}\rho^{\rm
III,3}_{g}\big{)}^{i_{3}}\big{(}\tilde{\rho}^{\rm
III,4}_{1}\big{)}^{i_{4}}\big{(}\tilde{\rho}^{\rm
III,5}_{1}\big{)}^{i_{5}}\big{(}\tilde{\rho}^{\rm
III,6}_{2}\big{)}^{i_{6}}\big{(}\tilde{\rho}^{\rm
III,7}_{2}\big{)}^{i_{7}}}\,,$ (170)

where the propagators are defined by

$\displaystyle\rho^{\rm III,1}_{g}$ $\displaystyle=\ell_{1}^{2}\,,$
$\displaystyle\quad\rho^{\rm III,2}_{g}$ $\displaystyle=\ell_{2}^{2}\,,$
$\displaystyle\quad\rho^{\rm III,3}_{g}$
$\displaystyle=(\ell_{1}+\ell_{2}+q)^{2}\,,$ (171)
$\displaystyle\tilde{\rho}^{\rm III,4}_{1}$ $\displaystyle=2\ell_{1}\cdot
u_{1}\,,$ $\displaystyle\quad\tilde{\rho}^{\rm III,5}_{1}$
$\displaystyle=-2\ell_{2}\cdot u_{1}\,,$ $\displaystyle\quad\tilde{\rho}^{\rm
III,6}_{2}$ $\displaystyle=-2\ell_{1}\cdot u_{2}\,,$
$\displaystyle\tilde{\rho}^{\rm III,7}_{2}$ $\displaystyle=2\ell_{2}\cdot
u_{2}\,,$ $\displaystyle\qquad\rho^{\rm III,8}_{g}$
$\displaystyle=(\ell_{1}+q)^{2},$ $\displaystyle\qquad\rho^{\rm III,9}_{g}$
$\displaystyle=(\ell_{2}+q)^{2}\;.$

After the IBP reduction, we have ten master integrals defined by

$\displaystyle f^{\text{III}}_{1}$
$\displaystyle=-q^{2}\epsilon^{2}j^{\text{III}}_{2,2,1,0,0,0,0,0,0}\,,$ (172)
$\displaystyle f^{\text{III}}_{2}$
$\displaystyle=-\sqrt{y^{2}-1}\epsilon^{4}j^{\text{III}}_{1,1,1,1,0,0,1,0,0}\,,$
$\displaystyle f^{\text{III}}_{3}$
$\displaystyle=q^{2}\sqrt{y^{2}-1}\epsilon^{3}j^{\text{III}}_{1,1,2,1,0,0,1,0,0}\,,$
$\displaystyle f^{\text{III}}_{4}$
$\displaystyle=q^{2}y\epsilon^{3}j^{\text{III}}_{1,1,2,1,0,0,1,0,0}+q^{2}\epsilon^{2}j^{\text{III}}_{1,1,1,2,0,0,2,0,0}\,,$
$\displaystyle f^{\text{III}}_{5}$
$\displaystyle=q^{2}\sqrt{y^{2}-1}\epsilon^{3}j^{\text{III}}_{1,1,2,0,1,0,1,0,0}\,,$
$\displaystyle f^{\text{III}}_{6}$
$\displaystyle=-(1-6\epsilon)\epsilon^{3}j^{\text{III}}_{1,1,1,1,1,0,0,0,0}\,,$
$\displaystyle f^{\text{III}}_{7}$
$\displaystyle=-q^{2}\left(y^{2}-1\right)\epsilon^{4}j^{\text{III}}_{1,1,1,1,1,1,1,0,0}\,,$
$\displaystyle f^{\text{III}}_{8}$
$\displaystyle=-|q|\epsilon^{3}j^{\text{III}}_{1,1,2,0,1,0,0,0,0}\,,$
$\displaystyle f^{\text{III}}_{9}$
$\displaystyle=|q|\epsilon^{3}j^{\text{III}}_{1,1,1,1,0,0,2,0,0}\,,$
$\displaystyle f^{\text{III}}_{10}$
$\displaystyle=|q|\sqrt{y^{2}-1}\epsilon^{4}j^{\text{III}}_{1,1,1,1,1,0,1,0,0}\,.$

Differentiating $f^{\rm III}_{i}$ with respect to $y$ and changing the
variable to $x$, one can derive a set of differential equations in the
canonical form

$\frac{\mathrm{d}\vec{f}_{\rm III}}{\mathrm{d}x}=\epsilon\bigg{(}\frac{A_{\rm
III,0}}{x}+\frac{A_{\rm III,1}}{x-1}+\frac{A_{\rm
III,-1}}{x+1}\bigg{)}\vec{f}_{\rm III}\,,\qquad\vec{f}_{\rm
III}=\big{\\{}f^{\rm III}_{1},f^{\rm III}_{2},\cdots,f^{\rm III}_{10}\\}\,,$
(173)

where $A_{\rm III}$ are $10\times 10$ constant matrices, which are divided
into two parts depending on even and odd powers in $|q|$

$A_{\mathrm{III},i}=\left(\begin{array}[]{cc}A_{\mathrm{III},i}^{(\mathrm{e})}&0\\\
0&A_{\mathrm{III},i}^{(\mathrm{o})}\end{array}\right)\,,$ (174)

and

$\displaystyle
A_{\mathrm{III},0}^{(\mathrm{e})}=\left(\begin{array}[]{ccccccc}0&0&0&0&0&0&0\\\
-\frac{1}{2}&-6&0&-1&0&0&0\\\ -\frac{3}{2}&0&2&-2&0&0&0\\\ 0&12&2&0&0&0&0\\\
-\frac{3}{4}&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0\\\
0&0&1&0&-2&0&0\end{array}\right)\,,\qquad A_{\mathrm{III},\pm
1}^{(\mathrm{e})}=\begin{pmatrix}0&0&0&0&0&0&0\\\ 0&6&0&0&0&0&0\\\
0&0&-2&0&0&0&0\\\ 0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0\\\ 0&0&0&0&0&0&0\\\
0&0&0&0&0&0&0\end{pmatrix}\,,$ (175)

and

$A_{\mathrm{III},0}^{(\mathrm{o})}=\left(\begin{array}[]{ccc}0&0&0\\\
0&-2&0\\\ 0&1&0\end{array}\right),\quad
A_{\mathrm{II},+1}^{(\mathrm{o})}=\left(\begin{array}[]{ccc}0&0&0\\\ 3&6&0\\\
0&0&0\end{array}\right),\quad
A_{\mathrm{III},-1}^{(\mathrm{o})}=\left(\begin{array}[]{ccc}0&0&0\\\
-3&-2&0\\\ 0&0&0\end{array}\right)\,.$ (176)

We present the solution of the differential equation under our convention in
(291).

We perform IBP reduction to the loop integrands for the $\rm III$ sector,
which are obtained by solving the off-shell recursions after the soft
expansion. Consequently, the amplitude $\mathcal{M}^{\rm III}$ is expanded in
terms of the master integrals (172)

$\mathcal{M}^{\rm III}=\sum_{i=1}^{10}c^{\rm III}_{i}f^{\rm
III}_{i}\Big{|}_{\text{non-analytic terms}}\,.$ (177)

The explicit form of the coefficients $c^{\rm III}_{i}$ and the value of the
master integrals $f^{\rm III}_{i}$ are listed in Appendix D.1.1 and (291)
respectively.

The classical amplitude $\mathcal{M}^{\rm III}$ can be obtained as

$\displaystyle\mathcal{M}^{\rm III}$
$\displaystyle=\frac{2\pi^{2}m_{1}^{4}m_{2}^{4}(2\sigma^{2}-1)^{3}}{3q^{2}(\sigma^{2}-1)}\Bigg{[}\pi^{2}-4\log^{2}x-\frac{6}{\epsilon^{2}}\Bigg{]}$
(178) $\displaystyle\
+\frac{4i\pi^{3}m_{1}^{3}m_{2}^{3}(m_{1}{+}m_{2})(2\sigma^{2}-1)}{|q|\sqrt{\sigma^{2}-1}}\Bigg{[}(4\sigma^{3}-2\sigma+1)\bigg{(}\frac{1}{\epsilon}-2\log
2\bigg{)}-\frac{(1-2\sigma^{2})^{2}}{\sigma-1}\Bigg{]}$ $\displaystyle\
+\frac{\pi^{2}m_{1}^{2}m_{2}^{2}}{15\epsilon}\Bigg{[}10m_{1}m_{2}\left(\frac{72\sigma^{7}-112\sigma^{5}+44\sigma^{3}-\sigma}{(\sigma^{2}-1)^{2}}-\frac{8\left(6\sigma^{4}-5\sigma^{2}+1\right)\log
x}{\sqrt{\sigma^{2}-1}}\right)$
$\displaystyle\qquad\qquad\quad+\left(m_{1}^{2}{+}m_{2}^{2}\right)\left(\frac{520\sigma^{8}{-}1272\sigma^{6}{+}1216\sigma^{4}{-}516\sigma^{2}{+}67}{(\sigma^{2}-1)^{2}}-\frac{80\sigma\log
x}{\sqrt{\sigma^{2}-1}}\right)\Bigg{]}\,.$

Next, we consider the $\overline{\rm III}$, the crossed counterpart of the
III-sector. The integral family and its master integrals of the $\overline{\rm
III}$ sector can be derived by the crossing relation from the III sector in
(182) and (172). The differential equations are also modified by the same
crossing relation. The explicit form of the master integrals is presented in
(292).

Using IBP reduction, we may represent the classical amplitude
$\mathcal{M}^{\overline{\rm III}}$ as a linear combination of the master
integrals $f^{\overline{\rm III}}_{i}$

$\mathcal{M}^{\overline{\rm III}}=\sum_{i=1}^{10}c^{\overline{\rm
III}}_{i}f^{\overline{\rm III}}_{i}\Big{|}_{\text{non-analytic terms}}\,.$
(179)

Similar in (177), we keep the non-analytic terms in $q$ or $\epsilon$, which
are relevant in the classical long-range interaction. The coefficients
$c^{\overline{\rm III}}_{i}$ are denoted in Appendix D.1.2. Note that the
coefficients for the III and $\overline{\rm III}$ sectors are in the crossing
relation,

$c^{\rm III}_{i}=c^{\overline{\rm
III}}_{i}\Big{|}_{\eqref{eq:crossed_diagram}}\,.$ (180)

On the other hand, we can derive the coefficients from the IBP reduction of
the loop integrands that are the solution of the off-shell recursions instead
of the crossing relation. This relation provides a consistency check.

Substituting the values of the coefficients and the master integrals into
(179), we have