# Two-photon exchange in leptophilic dark matter scenarios

Raghuveer Garani<EMAIL_ADDRESS>INFN Sezione di Firenze, Via G. Sansone 1,
I-50019 Sesto Fiorentino, Italy  Federico Gasparotto
<EMAIL_ADDRESS>Dipartimento di Fisica e Astronomia, Università
di Padova, Via Marzolo 8, 35131 Padova, Italy  INFN, Sezione di Padova, Via
Marzolo 8, 35131 Padova, Italy  Pierpaolo Mastrolia
<EMAIL_ADDRESS>Dipartimento di Fisica e Astronomia, Università
di Padova, Via Marzolo 8, 35131 Padova, Italy  INFN, Sezione di Padova, Via
Marzolo 8, 35131 Padova, Italy
Henrik J. Munch<EMAIL_ADDRESS>Dipartimento di Fisica e Astronomia,
Università di Padova, Via Marzolo 8, 35131 Padova, Italy  INFN, Sezione di
Padova, Via Marzolo 8, 35131 Padova, Italy  Sergio Palomares-Ruiz
<EMAIL_ADDRESS>Instituto de Fìsica Corpuscular (IFIC), CSIC-Universitat
de València, Apartado de Correos 22085, E-46071 València, Spain  Amedeo Primo
<EMAIL_ADDRESS>Department of Physics, University of Zürich, CH-8057
Zürich, Switzerland

###### Abstract

In leptophilic scenarios, dark matter interactions with nuclei, relevant for
direct detection experiments and for the capture by celestial objects, could
only occur via loop-induced processes. If the mediator is a scalar or pseudo-
scalar particle, which only couples to leptons, the dominant contribution to
dark matter-nucleus scattering would take place via two-photon exchange with a
lepton triangle loop. The corresponding diagrams have been estimated in the
literature under different approximations. Here, we present new analytical
calculations for one-body two-loop and two-body one-loop interactions. The
two-loop form factors are presented in closed analytical form in terms of
generalized polylogarithms up to weight four. In both cases, we consider the
exact dependence on all the involved scales, and study the dependence on the
momentum transfer. We show that some previous approximations fail to correctly
predict the scattering cross section by several orders of magnitude. Moreover,
we show that form factors, in the range of momentum transfer relevant for
local galactic dark matter, are smaller than their value at zero momentum
transfer, which is usually considered.

††preprint: ZU-TH 21/21

##  I Introduction

Precise perturbative calculations in quantum field theory, within the Standard
Model (SM) and in theories beyond the SM (BSM), require the computation of
loop integrals. Over the past decades, many elaborate computational tools for
the calculation of multi-loop Feynman diagrams have been developed (see for
reviews, e.g., Refs. [1, 2, 3, 4, 5]) and have been applied for the
calculation of many processes.

The diagrams in Fig. 1, with generic internal fermion and boson masses, often
contribute to SM and BSM observables. Within perturbative quantum
chromodynamics (QCD), the two-loop vertex graphs appear in second-order
corrections to the form factors involving heavy quarks, and have been
considered in the calculation of vector and axial-vector form factors in the
heavy quark forward-backward asymmetry [6, 7, 8, 9, 10, 11, 12] or of the
scalar and pseudo-scalar form factors in Higgs decays to heavy quarks [13, 14,
12, 15, 16, 17, 18]. Other static quantities, such as anomalous magnetic and
electric moments, also receive contributions from similar two-loop diagrams
[19, 20, 21, 22, 23, 24, 25] via, e.g., Barr-Zee diagrams [26] (involving
bosons with different masses in the loop). On the other hand, non-relativistic
QCD [27], which describes the dynamics of heavy quark-antiquark pairs at
energy scales well below their masses, provides accurate predictions for,
e.g., the decay rates and binding energies of heavy quarkonia and for top-pair
production near threshold in electron-positron annihilation [28, 29]. Two-loop
diagrams with the generic shape of those in Fig. 1 also appear in the
calculation of several matching coefficients between QCD and non-relativistic
QCD at two-loop order [30]. Within the context of the quark-loop model [31]
(known to be equivalent to vector meson dominance models [32]), these diagrams
represent the dominant contribution in the calculation of rare pseudo-scalar
meson decays into leptons [31, 33, 34].

Moreover, the decay of a scalar or pseudo-scalar boson into quarks has been
investigated in great detail, within the context of the SM Higgs boson or of
BSM non-standard Higgs bosons. Early computations of second-order corrections
to these decays have been obtained for the case of massless quarks [35] and in
expansions in terms of the ratio of the external quarks and boson masses [36,
37, 38, 39, 40, 41, 42]. Analytical expressions for the two-loop diagrams in
Fig. 1, where two photons are exchanged, were earlier derived for the cases
when the mass of the external fermion is either vanishing or equal to that of
the particles running in the triangle loop (equal-mass case) [43, 44, 45, 46,
47, 13, 14, 12]. More recently, analytical results for the different-mass case
have been obtained in the context of Higgs decays into bottom [16] and charm
quarks [18].

Figure 1: Two-boson exchange diagrams for interactions between a scalar or
pseudo-scalar boson and two external fermions,
$\phi(q)\,\bar{\psi}(p^{\prime})\,\Lambda_{S/P}\,\psi(p)$
($\Lambda_{S}=\mathbb{I}$ and $\Lambda_{P}=\gamma_{5}$), via a lepton triangle
loop which exchanges two bosons with an external fermion: one-body two-loop
interactions and two-body one-loop interactions.

Interestingly, the same graphs are also encountered in theoretical scenarios
of particle dark matter (DM). Indeed, the diagrams shown in Fig. 1 can
represent the dominant contribution for interactions with nucleons of a
particular class of DM candidates, which have no tree-level couplings to
quarks, but only to leptons. These leptophilic scenarios have been proposed to
reduce the tension of the DM interpretation among different direct detection
experiments [48, 49, 50, 51, 52, 53, 54, 55], to solve some cosmic-ray
anomalies [56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67], to explain the
gamma-ray galactic center excess [68], to alleviate the tension in the
measured muon magnetic moment [69, 70, 71], as potential signals in collider
searches [72, 73, 74, 75, 76, 77] or in gravitational wave detectors [75], as
a way to capture DM in celestial objects [50, 78, 79, 80, 81, 82, 83, 84] and,
within the context of sub-GeV DM, to discuss future strategies of direct
detection [85, 86, 87, 88], as a source of distortion of the cosmic-ray
electron spectrum [89] or as a way to produce DM up-scatterings to
relativistic energies [90, 91, 92].

If the DM particle is leptophilic, its tree-level coupling to fermions
involves just leptons, but not quarks. Therefore, the scattering cross section
off nuclei (or off quarks) could only proceed via loop-induced interactions.
In the case of a scalar or pseudo-scalar mediator, the one-body one-loop
contribution (with a lepton loop) to the DM-quark/nucleon coupling is
suppressed by the masses of electroweak gauge bosons and the leading
interactions with nuclei arise from diagrams as those in Fig. 1, with two
photons exchanged with external protons. The two left diagrams represent one-
body two-loop interactions and the two right diagrams, two-body one-loop
processes. In such scenarios, a mediator particle connects DM to SM particles,
so DM-nucleon elastic scattering proceeds via a lepton triangle loop and two
photons. Thus, in some leptophilic scenarios, those diagrams could govern the
capture of DM in celestial objects and the detection rates in terrestrial
experiments.

There have been attempts to compute the amplitude of the above diagrams in the
context of detection of leptophilic DM. In Refs. [50, 93], an effective
operator product expansion approach was discussed. After integrating out the
lepton loop, an effective interaction term $\bar{\chi}\chi
F_{\mu\nu}F^{\mu\nu}$ was obtained (with $\chi$ being the DM field).111This
coupling is analogous to that in Rayleigh DM scenarios [94], although
different UV completions can give rise to it without a fermion triangle loop
[95, 96]. In Ref. [50], this was then matched onto the external nucleus
current to obtain an effective DM-nucleus interaction in the heavy-nucleon
limit, which actually represents a two-body process [97]. However, this
approximation is not reliable when the momentum transfer is comparable to the
fermion masses in the triangle loop [50]. Furthermore, we show that this is
not appropriate in general. On the other hand, in Ref. [80], direct
integration of the one-body two-loop amplitude is performed in the limit of
vanishing momentum transfer and an approximate analytical expression is
reported.

In this work, we present the full analytical calculation of the one-body two-
loop and two-body one loop diagrams (Fig. 1), keeping the complete dependence
on the kinematic variables, namely the momentum transfer of the mediator and
the two fermion masses. We consider Feynman diagrams in dimensional
regularization, with $d=4-2\epsilon$ dimensions (for a small, positive
$\epsilon$ parameter), and exploit the integration-by-parts identities (IBPs)
[98], to express the scalar and pseudo-scalar form factors in terms of a set
of independent master integrals (MIs) [99, 100], recently computed in [16,
18]. Then, we consider the four dimensional limit, by taking $\epsilon\to 0$,
where the form factors are finite, and express the results of the one-body
diagrams in terms of generalized polylogarithms (GPLs) [101, 102, 103, 104,
105, 106]. Besides the expressions corresponding to general kinematics, we
provide the analytic expressions of the form factors in two limits: the soft
limit, for vanishing momentum transfer of the mediator; and the equal-mass
limit, for internal fermions with equal masses. We use the novel results of
the two-loop one-body and of the one-loop two-body diagrams to investigate the
leading contribution to the scattering of leptophilic DM off a few sample of
nuclei, like helium, oxygen, iron and xenon, relevant for direct DM detection
experiments and for the capture of DM by celestial objects. Finally, we
critically compare our results with previous estimates, finding that the exact
calculation of the form factors presented here turns out to be crucial for
phenomenological analyses, and that naive approximations of the loop integrals
might yield incorrect conclusions.

This article is organized as follows. We derive the relevant form factors for
scalar and pseudo-scalar interactions for the one-body and the two-body
processes, write down the amplitudes and obtain the analytical solution for
the form factors of one-body interactions in terms of GPLs in Section II.1.
The calculation of the contribution to the amplitude from two-body
interactions is described in Section II.2. Thus, Section II represents the
core of this paper, where we describe and discuss the main results. We compare
our results to the approximate results available in the literature in Section
III. Finally, in Section IV we draw our conclusions. Details about the
calculations of two-loop diagrams are presented in Appendices A, B and C.

##  II Amplitudes and form factors

We consider the case of a DM particle ($\chi$) that couples to a SM singlet
scalar ($\phi_{S}$) or pseudo-scalar ($\phi_{P}$), which in turn, couples to
SM leptons ($\ell$) at tree level, but not to quarks, through interaction
terms of the form

$-{\cal L}_{S}\supset
g_{S}\,\phi_{S}\,\bar{\ell}\,\ell+g_{\chi}\,\phi_{S}\,\Gamma_{\chi}\hskip
28.45274pt\textrm{or}\hskip 28.45274pt-{\cal L}_{P}\supset
i\,g_{P}\,\phi_{P}\,\bar{\ell}\,\gamma_{5}\,\ell+g_{\chi}\,\phi_{P}\,\Gamma_{\chi}~{}.$
(1)

where $\Gamma_{\chi}=\\{\chi^{\dagger}\chi\\}$ (scalar DM) or
$\\{\bar{\chi}\,\chi,\bar{\chi}\,\gamma_{5}\,\chi\\}$ (fermion DM). The
leading interactions with nucleons are induced by the exchange of two virtual
photons. In general, the two virtual photons can scatter off a single nucleon
[50, 93] or two nucleons [94, 97], and these two types of processes govern the
elastic scattering of a DM particle off a nucleus $A$, with the total
amplitude given by the sum of their individual amplitudes,

${\cal M}_{\beta}=\langle A(p^{\prime})|\,\mathcal{L}\,|A(p)\rangle={\cal
M}^{(1)}_{\beta}+{\cal M}^{(2)}_{\beta}~{},$ (2)

with $\beta=\\{S,P\\}$ for scalar or pseudo-scalar interactions.

The computation of these contributions is the main goal of this paper and we
discuss them next.

###  II.1 One-body interactions

The leading one-body interactions with quarks/nucleons are induced at two-loop
level through the two left diagrams shown in Fig. 1, where SM charged leptons
that run in the triangle loop exchange two photons with an external
quark/nucleon. As mentioned above, our main target is the scattering of non-
relativistic DM particles off nucleons and nuclei. In principle, one could
proceed in terms of interactions with quarks and then match to the effective
theory that would describe the interaction of non-relativistic DM with non-
relativistic nucleons. However, in order to do so, the dressed quark
propagators (and vertices) within the nucleon would have to be used. Thus, for
the low energies of interest, for which the internal structure of nucleons is
not probed, it is more convenient to directly consider the well-known coupling
of photons to nucleons. This is given by the electromagnetic current for a
nucleon, which is generically written as

$\langle
N(p^{\prime})|J^{\mu}(q^{2})|N(p)\rangle=\bar{u}_{N}(p^{\prime})\,\left(\gamma^{\mu}\,F_{1}(q^{2})+i\,\sigma^{\mu\nu}\frac{q_{\nu}}{2m_{N}}\,F_{2}(q^{2})\right)\,u_{N}(p)~{},$
(3)

where $q^{2}=(p^{\prime}-p)^{2}\equiv t$ is the squared of the exchanged
momentum, $m_{N}$ is the nucleon mass, and $F_{1}(q^{2})$ and $F_{2}(q^{2})$
are the Dirac and Pauli form factors. At $q^{2}=0$, $F_{1}(0)$ represents the
electric charge of the nucleon and $F_{2}(0)$ its anomalous magnetic moment.
In the non-relativistic limit (i.e., $q/m_{N}<1$) relevant for the problem
under discussion, $F_{1}(q^{2})\simeq F_{1}(0)$ and $F_{2}(q^{2})\simeq
F_{2}(0)$. Moreover, the term involving $F_{2}(q^{2})$ is proportional to
$q/m_{N}$ and thus, can also be neglected at the lowest order. Therefore, in
the non-relativistic limit, the electromagnetic coupling of photons to
nucleons is analogous to that of quarks, the only difference being the charge.
Obviously, at this order, only protons contribute. For the sake of comparison
with other works, however, we will make this simplification for the photon-
nucleon vertex, but we will keep the rest of the treatment fully relativistic.

The amplitude of the two-loop vertex diagrams in Fig. 1 can be written, in
terms of the interaction with nucleons, as

$\tilde{\mathcal{M}}_{\beta}^{(1)}(t;m_{N})=i\,g_{\beta}\,Q_{N}^{2}\,\sum_{\ell=e,\mu,\tau}Q_{\ell}^{2}\,\left(\bar{u}_{N}(p^{\prime})\,\Gamma_{\beta}(t;m_{N},m_{\ell})\,u_{N}(p)\right)~{},$
(4)

where $m_{\ell}$ ($Q_{\ell}=1$) and $m_{N}$ ($Q_{N}=\\{1,0\\}$, for $N={p,n}$)
are the masses (electric charges) of lepton $\ell$ and nucleon $N$,
respectively, and $p^{2}=p^{\prime 2}=m_{N}^{2}$. The operator
$\Gamma_{\beta}(t,m_{N},m_{\ell})$ results from the sum of the two diagrams,
which give identical contributions,

$\Gamma_{\beta}(t;m_{N},m_{\ell})=-32\,\pi^{2}\,\alpha_{\rm
em}^{2}\,\int\frac{{\differential^{d}k_{1}}}{(2\pi)^{4}}\frac{{\differential^{d}k_{2}}}{(2\pi)^{4}}\frac{g_{\mu\rho}\,\gamma^{\rho}\,\left(\not{k_{2}}-m_{N}\right)\,g_{\nu\sigma}\,\gamma^{\sigma}\,\text{Tr}^{\mu\nu}_{\beta}(q,k_{1},k_{2};m_{N},m_{\ell})}{D_{1}\,D_{2}\,D_{3}\,D_{4}\,D_{5}\,D_{6}}~{},$
(5)

where $\alpha_{\rm em}=e^{2}/(4\,\pi)$ and $\text{Tr}^{\mu\nu}_{\beta}$ is the
Dirac trace

$\text{Tr}^{\mu\nu}_{\beta}(t,k_{1},k_{2};m_{N},m_{\ell})=\text{Tr}\Big{\\{}\left(\not{k_{1}}+\not{q}+m_{\ell}\right)\Lambda_{\beta}\left(\not{k_{1}}+m_{\ell}\right)\gamma^{\mu}\left(\not{k_{1}}+\not{k_{2}}+\not{p}^{\prime}+m_{\ell}\right)\gamma^{\nu}\Big{\\}}~{},$
(6)

with $\Lambda_{S}=\mathbb{I}$ and $\Lambda_{P}=\gamma_{5}$. The inverse
propagators $D_{i}$ in Eq. (5) are

$\displaystyle D_{1}=k_{1}^{2}-m_{\ell}^{2}~{},\quad
D_{2}=(k_{1}+q)^{2}-m_{\ell}^{2}~{},\quad D_{3}=k_{2}^{2}-m_{N}^{2}~{},$
$\displaystyle D_{4}=(k_{2}+p)^{2}~{},\quad
D_{5}=(k_{2}+p^{\prime})^{2}~{},\quad
D_{6}=(k_{1}+k_{2}+p^{\prime})^{2}-m_{\ell}^{2}~{}.$ (7)

From Lorentz invariance, the operators $\Gamma_{\beta}$ must have the form

$\Gamma_{\beta}(t;m_{N},m_{\ell})=A_{\beta}(t;m_{N},m_{\ell})\,\Lambda_{\beta}+B_{\beta}(t;m_{N},m_{\ell})\,\left(\not{p}^{\prime}+\not{p}\right)\,\Lambda_{\beta}+C_{\beta}(t;m_{N},m_{\ell})\,\left(\not{p}^{\prime}-\not{p}\right)\,\Lambda_{\beta}~{},$
(8)

where the form factors $A_{\beta}$, $B_{\beta}$ and $C_{\beta}$ are scalar
functions of the square of the momentum transfer, $t$, as well as of the
fermion masses $m_{N}$ and $m_{\ell}$. In order to write the effective
interaction term $\phi_{\beta}\bar{N}\Lambda_{\beta}N$ in the Lagrangian, we
use the Dirac equation for the external fermions. That is, we make the
following identifications:

$\displaystyle\bar{u}_{N}(p^{\prime})\,\Lambda_{\beta}\,u_{N}(p)$
$\displaystyle\to\bar{N}(p^{\prime})\,\Lambda_{\beta}\,N(p)~{},$
$\displaystyle\bar{u}_{N}(p^{\prime})\,(\not{p}^{\prime}+\not{p})\,\left(1+\gamma_{5}\right)\,u_{N}(p)$
$\displaystyle\to 2\,m_{N}\,\bar{N}(p^{\prime})N(p)~{},$
$\displaystyle\bar{u}_{N}(p^{\prime})\,(\not{p}^{\prime}-\not{p})\,\left(1+\gamma_{5}\right)\,u_{N}(p)$
$\displaystyle\to 2\,m_{N}\,\bar{N}(p^{\prime})\,\gamma_{5}\,N(p)~{}.$ (9)

In this way, this dark mediator-nucleon interaction can be written as

$\mathcal{L}_{\beta}\supset
g_{\beta}\,Q_{N}^{2}\,\left[\sum_{\ell=e,\mu,\tau}Q_{\ell}^{2}\,{\cal
F}_{\beta}^{\rm
1b}(t;m_{N},m_{\ell})\right]\phi_{\beta}\,\bar{N}\,\Lambda_{\beta}\,N~{},$
(10)

where ${\cal F}_{S}^{\rm 1b}=A_{S}+2\,m_{N}\,B_{S}$ and ${\cal F}_{P}^{\rm
1b}=A_{P}+2\,m_{N}\,C_{P}$.

Since we are interested in the elastic scattering of DM particles off nuclei
in the non-relativistic limit, we can write the one-body contribution to the
DM-nucleus amplitude, for small momentum transfer, as

${\cal M}^{(1)}_{\beta}(\bm{q})=\langle
A(p^{\prime})|\mathcal{L}_{\beta}|A(p)\rangle=G\,g_{\beta}\,{\cal N}\,\left[
\sum_{\ell=e,\mu,\tau}{\cal F}_{\beta}^{\rm 1b}(t;m_{N},m_{\ell})\right]
\,F^{(1)}\left(\frac{\bm{q}^{2}}{Q_{0}^{2}}\right)~{},$ (11)

where $G$ includes the $\phi_{\beta}$ propagator and the DM-$\phi_{\beta}$
coupling $g_{\chi}$, ${\cal N}$ depends on the type of interaction, and
$F^{(1)}\left(-\bm{q}^{2}/Q_{0}^{2}\right)$ is the nuclear form factor, which
takes into account the coherence of the interaction. In the non-relativistic
limit, it depends on the square of the three-momentum transfer
$\bm{q}^{2}\simeq-t$, and throughout this work, it is approximated as

$F^{(1)}(\bm{q}^{2}/Q_{0}^{2})=e^{-\bm{q}^{2}/Q_{0}^{2}}.$ (12)

The latter expression of the form factor can be adapted both to the scalar
(spin-independent) and to pseudo-scalar (spin-dependent) cases, by choosing a
proper definition of $Q_{0}$ (which would be different for the two cases). For
the scalar case, $1/Q_{0}$ represents the radius of the nucleus,
$Q_{0}=0.48\,(0.3+0.89\,A^{1/3})^{-1}$ GeV. Also note that, in the non-
relativistic limit, the nuclear amplitude only depends on the number of
protons in the nucleus, ${\cal N}=Z$, with $Z$ the atomic number of the
nucleus.

The form factors ${\cal F}_{\beta}^{\rm 1b}$ can be extracted from
$\Gamma_{\beta}$ by projection,

${\cal F}_{\beta}^{\rm 1b}(t;m_{N},m_{\ell})=\frac{1}{2\left(p^{\prime}\pm
p\right)^{2}}\,\text{Tr}\left\\{\Lambda_{\beta}\left(\not{p}\pm
m_{N}\right)\,\Gamma_{\beta}\,\left(\not{p}^{\prime}\pm
m_{N}\right)\right\\},$ (13)

where the sign $+$ corresponds to the scalar case ($\Lambda_{S}=\mathbb{I}$)
and the sign $-$ to the pseudo-scalar case ($\Lambda_{P}=\gamma_{5}$). The two
form factors are both UV and IR finite in four dimensions, so that Eq. (13)
could be directly evaluated in $d=4$. However, as we describe in the
following, it is technically more convenient to exploit the algebraic
properties of Feynman integral in dimensional regularization, and to recover
the four dimensional result at the end of the calculation. Fermionic traces
and loop integration are performed, therefore, in arbitrary $d$ dimensions, by
adopting the prescription of Refs. [107, 108, 109] for the $d$-dimensional
treatment of $\gamma_{5}$.

Upon the calculation of the fermion traces, the form factors ${\cal
F}_{\beta}^{\rm 1b}(t;m_{N},m_{\ell})$ are expressed as linear combinations of
two-loop Feynman integrals of the type

${\cal F}_{\beta}^{\rm
1b}(t;m_{N},m_{\ell})=\sum_{\vec{n}}c_{\vec{n},\beta}(m_{i}^{2},t,\epsilon)\int\frac{{\differential^{d}k_{1}}}{(2\pi)^{d}}\frac{{\differential^{d}k_{2}}}{(2\pi)^{d}}\,\frac{D_{7}^{n_{7}}}{D_{1}^{n_{1}}D_{2}^{n_{2}}D_{3}^{n_{3}}D_{4}^{n_{4}}D_{5}^{n_{5}}D_{6}^{n_{6}}}\,,\quad
n_{i}\in\mathbb{N}~{},$ (14)

where $c_{\vec{n}}$ are coefficients that depend on the kinematic variables
and on the dimensional regulator $\epsilon=(4-d)/2$. The auxiliary inverse
propagator $D_{7}=(k_{1}-p)^{2}$ has been introduced to complement the six
denominators defined in Eq. (7), so that all the seven possible scalar
products between the loop momenta $k_{1,2}$ and the independent external
momenta $q\,,p$ are expressible as linear combinations of $D_{i}$.

The exact expression of ${\cal F}_{S}^{\rm 1b}$ contains 34 distinct Feynman
integrals of the type in Eq. (14), that of ${\cal F}_{P}^{\rm 1b}$ only a
subset of 5 integrals. However, these 34 integrals are not all independent
from each other, due to linear relations entailed by $d$-dimensional IBPs. We
employ the packages Reduze2 [100] and LiteRed [110] for generating the IBPs
needed to reduce the form factors to a combination of a minimal set of 20
independent MIs.

The set of 20 MIs have been analytically computed through the differential
equation method in Ref. [16] in the context of the evaluation of top Yukawa
corrections to the decay of a Higgs boson into a $b\bar{b}$-pair, which
involve two-loop triangle diagrams similar to the one considered here.
Therefore, we adopt the same method of Ref. [16] to evaluate the MIs, as we
summarize in the following.

The integrals under consideration depend on three kinematic invariants, the
square of the momentum transfer $t$, and the nucleon and lepton masses,
$m_{N}$ and $m_{\ell}$. For a given $m_{\ell}$, the other two variables can be
combined into two dimensionless parameters $x$ and $y$, which retain, up to a
trivial scaling factor depending on $m_{\ell}$, the full dependence of the
integrals on the kinematics. For the problem under consideration, it is
convenient to define

$t=-m_{\ell}^{2}\,\frac{(1-x^{2})^{2}}{x^{2}}~{},\qquad
m_{N}^{2}=m_{\ell}^{2}\,\frac{(1-x^{2})^{2}}{(1-y^{2})^{2}}\frac{y^{2}}{x^{2}}~{}.$
(15)

Owing to IBPs, MIs fulfill systems of coupled first order differential
equations in $t/m_{\ell}^{2}$ and $m_{N}^{2}/m_{\ell}^{2}$, or equivalently,
in $x$ and $y$. We determine the analytic expression of the MIs by solving, as
a series expansion around $\epsilon=0$, the corresponding system of
differential equations in $x$ and $y$. The solution of such differential
equations can be greatly simplified with a suitable choice of the basis of
MIs. We employ the Magnus exponential algorithm [111, 112] to identify a basis
that obeys canonical differential equations in the sense of Ref. [113],

$\displaystyle\partial_{x}\mathbf{I}(\epsilon,x,y)$ $\displaystyle=$
$\displaystyle\epsilon\,\mathbb{A}_{x}(x,y)\,\mathbf{I}(\epsilon,x,y)~{},$
$\displaystyle\partial_{y}\mathbf{I}(\epsilon,x,y)$ $\displaystyle=$
$\displaystyle\epsilon\,\mathbb{A}_{y}(x,y)\,\mathbf{I}(\epsilon,x,y)~{},$
(16)

where $\mathbb{A}_{x}$ and $\mathbb{A}_{y}$ are rational matrices and
$\mathbf{I}_{\beta}(\epsilon,x,y)$ is a vector that collects the basis
integrals. The factorized dependence of the differential equations on the
dimensional regulator $\epsilon$ streamlines the derivation of their general
solutions in terms of well-known special functions, GPLs [101, 102, 103, 104,
105]. With the general solution of Eq. (II.1) at hand, a proper choice of the
boundary conditions, which is dictated by physical requirements, allows the
complete determination of the required MIs and, hence, provides all the
information information required to obtain the analytic expression of the form
factors. For the technical details of the computation of the MIs, we refer the
reader to Appendix A.

After combining the analytic expression of the MIs into the form factors, we
observe the expected analytic cancellation of all $\epsilon$-poles possessed
by the individual integrals and obtain a finite result for both form factors
in the four-dimensional limit $\epsilon\to 0$. The results for the scalar
(${\cal F}_{S}^{\rm 1b}$) and pseudo-scalar (${\cal F}_{P}^{\rm 1b}$) form
factors are presented in Fig. 2, as a function of the squared momentum
transfer $t$. We show the results for target protons, which, for the values of
the considered momenta, are the relevant degrees of freedom in the scattering
of local galactic DM off standard matter, and present the partial
contributions from the diagrams with the three charged leptons in the triangle
loop: $e$ (black dashed curves), $\mu$ (green dashed curves) and $\tau$ (blue
dashed curves); and the sum of the three contributions (red solid curves).
Note that, in the non-relativistic limit, $-t\simeq|\bm{q}|^{2}$, where
$\bm{q}$ is the three-momentum transfer. In this limit, the maximum momentum
transfer to a nucleus is $|\bm{q}|_{\rm max}=2\,\mu\,v_{\rm
rel}<m_{A}\,v_{0}\sim 100$ MeV, where $\mu$ is the reduced mass of the DM-
nucleus system, $m_{A}\sim 100$ GeV is a typical target mass in direct
detection experiments,222In the process of DM capture in the Sun or Earth, the
target nuclei are lighter, and so is the momentum transfer. For DM capture in
neutron stars, the heaviest target particles are neutrons, but the maximum
momentum transfer is $\mathcal{O}(500)$ MeV. $v_{\rm rel}$ is the relative
velocity and $v_{0}\sim 10^{-3}$ is the typical velocity of the DM halo in the
galactic vicinity. For larger values, the velocity distribution is expected to
fall off exponentially, so momentum transfers larger than $|\bm{q}|_{\rm max}$
are unlikely.

Figure 2: Form factors for one-body interactions of DM with protons, as a
function of the square of the exchanged momentum, $-t$: scalar form factor,
${\cal F}_{S}^{\rm 1b}(t;m_{p},m_{\ell})$ (left panel) and pseudo-scalar form
factor, ${\cal F}_{P}^{\rm 1b}(t;m_{p},m_{\ell})$ (right panel). We show the
results for interactions with a given lepton in the triangle loop: $e$ (black
dashed curves), $\mu$ (green dashed curves) and $\tau$ (blue dashed curves);
and for the total contribution (red solid curves). The gray shaded area
indicates the range of momentum transfer relevant for scattering of local
galactic DM off standard matter.

Let us first discuss the behavior of the form factors at vanishing square
momentum transfer ($t=0$), which is the usual limit adopted when studying
local DM scattering, in the calculation of nuclear recoils within direct
detection experiments, and in the calculation of capture of DM particles by
celestial objects. At zeroth order in the momentum transfer, the one-body
interaction only takes place with protons, so we will use $m_{N}=m_{p}$. For
the scalar form factor, we have obtained a closed expression, which reads

$\displaystyle{\cal F}_{S}^{\rm
1b}\left(t=0;z_{p\ell}\equiv\frac{m_{p}}{m_{\ell}}\right)$ $\displaystyle=$
$\displaystyle-\frac{2\,\alpha_{\rm
em}^{2}}{\pi^{2}\,z_{p\ell}}\Big{\\{}1-\frac{\ln(z_{p\ell})}{2}+f_{S}(z_{p\ell})+f_{S}(-z_{p\ell})\Big{\\}}~{},$
(17) $\displaystyle\textrm{with}\hskip 19.91692ptf(\tau)$
$\displaystyle\equiv$
$\displaystyle\frac{1}{4\,\tau^{2}}\left(4+3\,\tau+\tau^{3}\right)\Big{[}\ln|\tau|\ln(1+\tau)+{\rm
Li}_{2}(-\tau)\Big{]}~{}.$

Details of this calculation are provided in Appendix B.

Two interesting limits can be studied in terms of the ratio of the proton mass
and the mass of the leptons in the triangle loop, $z_{p\ell}\equiv
m_{p}/m_{\ell}$. For lepton masses larger than the proton mass, $z_{p\ell}\ll
1$,

${\cal F}_{S}^{\rm 1b}(t=0;z_{p\ell}\ll 1)\simeq\frac{\alpha_{\rm
em}^{2}}{\pi^{2}}\,z_{p\ell}\,\left[\,\ln\left(z_{p\ell}^{-1}\right)+\frac{13}{12}\right]~{}.$
(18)

This limit only (approximately) applies to the tau lepton loop. For the
electron and muon loops, one should consider the opposite limit, $z_{p\ell}\gg
1$, for which the scalar form factor approximates as

${\cal F}_{S}^{\rm 1b}(t=0;z_{p\ell}\gg 1)\simeq\frac{\alpha_{\rm
em}^{2}}{4}\,\left(1+\frac{3}{z_{p\ell}^{2}}\right)~{}.$ (19)

At first order, this limit is independent of the proton-lepton mass ratio and
provides an accurate estimate for the electron- and muon-loop contributions at
$t=0$. It also helps to understand the behavior of the contribution in the
opposite limit, Eq. (18), which is larger the larger the ratio $z_{p\ell}$,
but it has a limiting value given by Eq. (19). In short, it turns out that for
small $|t|$, the electron and muon loops contribute approximately the same to
the scalar form factor, whereas the tau-loop contribution is smaller, as can
be seen from Fig. 2.

For the pseudo-scalar form factor, an analytical expression for the soft
asymptotic limit is presented in Appendix C. The result varies little at small
$|t|$ in the region depicted in Fig. 2, and the contributions from the three
lepton loops are larger than in the scalar case. Nevertheless, unlike the
scalar case, there is no analogous saturation behavior in the soft limit, but
the pseudo-scalar form factor slowly diverges as $t\to 0$. In this limit,
$\cal{F}_{P}^{\rm 1b}$ scales as $z_{p\ell}\simeq m_{p}/m_{\ell}$ (see
Appendix C), so the muon- and tau-loop contributions are always smaller than
the electron one. This is in contrast to the results at larger $|t|$ for any
of the two form factors, as we shall see next.

Let us now consider the behavior of the form factors as a function of the
square of the exchanged momentum. In both, the scalar and pseudo-scalar cases,
the form factor starts decreasing above $|t|\gtrsim m_{\ell}^{2}$. As a
result, around $|t|\sim(10^{-4}-10^{-1})~{}\textrm{GeV}^{2}$, the relevant
range of exchanged momentum for scatterings of local galactic DM off
nucleons/nuclei (gray area), the contribution from the electron lepton loop is
suppressed with respect to the other two. In that range, the scalar form
factor is approximately constant and equal to the sum of the contributions
from the diagrams with the muon and tau loop at $t=0$, being the former the
larger one.

The behavior of the pseudo-scalar form factor is qualitatively similar to the
scalar case, although the former is one-two orders of magnitude larger than
the latter in the relevant range of $t$ for local galactic DM-nucleon
scattering. The contribution of each lepton loop decreases from a similar
value of $t$ as in the scalar case, although the dependence on the lepton mass
in the loop at small $|t|$ is different, and for the pseudo-scalar form
factor, the diagram with the electron loop is the largest one at $|t|\lesssim
10^{-4}~{}\textrm{GeV}^{2}$. As the muon-loop contribution is larger than the
tau-loop one up to $|t|\sim(1-10)~{}\textrm{GeV}^{2}$, the former is the most
important one in the relevant range of $t$ (gray area in Fig. 2). Moreover,
the muon-loop contribution does not change much up to $|t|\gtrsim
10^{-1}~{}\textrm{GeV}^{2}$, so its value is similar to that at small $|t|$.
Note, however, that at small $|t|$ the dominant contribution is always that
with the electron loop. Thus, similarly to the scalar case, using the value of
the form factor at small $|t|$ is not a good approximation for DM scattering
off nucleons/nuclei. Finally, note that the pseudo-scalar coupling would
induce spin-dependent interactions which, in the non-relativistic limit,
reduce to a one-nucleon potential $\propto\bm{q}\cdot\bm{s}$, and thus, depend
on the spin content of the nucleus.

Finally, let us note that we have also performed the loop integrals, Eq. (5),
using the SecDec code [114, 115, 116, 117] and we have checked that the
numerical results so obtained match very well with the analytical results
presented above, in regions where we observe numerical cancellation of all
$\epsilon$-poles. This occurs for relatively large momentum transfers, larger
than the the lepton mass (and for larger momentum transfers for the scalar
case). For relatively small momentum transfers, however, numerical results
from SecDec do not converge at the required precision, so we cannot reliably
compare them to the analytical results we obtain in those regions.

###  II.2 Two-body interactions

Let us now consider two-body DM interactions within a nucleus for the scalar
case. Within leptophilic scenarios, for local DM-nucleon scatterings, lepton
masses could be comparable to the typical momentum transfer. Thus, we need to
retain the full form of the loop function without integrating out heavy
leptons and matching onto the Rayleigh operator, $\bar{\chi}\chi
F_{\mu\nu}F^{\mu\nu}$. The amplitude of the two-body two-boson exchange
diagrams in Fig. 1 (including the external DM particles) takes the form

${\cal
M}_{S,ij}=\sum_{\ell=e,\mu,\tau}\frac{\Delta_{S}^{\mu\nu}(q_{i},q_{j},m_{\ell})}{q^{2}_{i}\,q^{2}_{j}}\,\left(\bar{u}_{i}\gamma_{\mu}u_{i}\,\bar{u}_{j}\gamma_{\nu}u_{j}\right)~{},$
(20)

where $u_{i}$ and $u_{j}$ are the two proton spinors, and
$\Delta_{S}^{\mu\nu}(q_{i},q_{j},m_{\ell})$ is the result of the loop
integrals, which are summed up. Note that $q_{i}$ above are 4-momenta of the
photons that connect the vertices of the loop to the external protons, so that
$q=q_{i}+q_{j}$. The function $\Delta_{S}^{\mu\nu}(q_{i},q_{j},m_{\ell})$ has
the generic form

$\Delta_{S}^{\mu\nu}(q_{i},q_{j},m_{\ell})=\alpha\,g^{\mu\nu}+\beta\,q^{\mu}_{i}\,q^{\nu}_{j}+\gamma\,q^{\mu}_{j}\,q^{\nu}_{i}+\delta\,q^{\mu}_{j}\,q^{\nu}_{j}+\kappa\,q^{\mu}_{i}\,q^{\nu}_{i}~{},$
(21)

where the coefficients $\alpha,\beta,\gamma,\delta,\kappa$ are functions of
the photon momenta $q_{i}$, the external momentum transfer $q$ and the lepton
masses $m_{\ell}$. In the non-relativistic limit, up to
$\mathcal{O}(\bm{q}^{2}/m_{p}^{2})$, the term proportional to $g^{00}$ is the
most important one,333We have checked numerically that the contributions from
the terms proportional to $\beta,\gamma,\delta,\kappa$ are, indeed, smaller.
and the amplitude then reduces to

${\cal
M}_{S,ij}=G\,g_{S}\,\sum_{\ell=e,\mu,\tau}\frac{\bm{q}_{i}\cdot\bm{q}_{j}}{\bm{q}_{i}^{2}\,\bm{q}_{j}^{2}}\,V_{\alpha}(\bm{q}_{i},\bm{q}_{j},m_{\ell})\,\left(1+{\cal
O}(\bm{q}^{2}/m_{p}^{2})\right)~{},$ (22)

where,

$\displaystyle V_{\alpha}(\bm{q}_{i},\bm{q}_{j},m_{\ell})$ $\displaystyle=$
$\displaystyle-
V_{\ell}\,\frac{3\,m_{\ell}^{2}}{\bm{q}_{i}\cdot\bm{q}_{j}}\left\\{1+\frac{\left(\bm{q}^{2}_{j}-\bm{q}^{2}_{i}+\bm{q}^{2}\right)}{\lambda(\bm{q}^{2}_{i},\bm{q}^{2}_{j},\bm{q}^{2})}\,\Omega(\bm{q}_{i}^{2},m_{\ell})+\frac{\left(\bm{q}^{2}_{i}-\bm{q}^{2}_{j}+\bm{q}^{2}\right)}{\lambda(\bm{q}^{2}_{i},\bm{q}^{2}_{j},\bm{q}^{2})}\,\Omega(\bm{q}_{j}^{2},m_{\ell})\right.$
$\displaystyle\hskip
62.59596pt+\left(\frac{1}{\bm{q}^{2}}+\frac{\bm{q}^{2}_{i}+\bm{q}^{2}_{j}-\bm{q}^{2}}{\lambda(\bm{q}^{2}_{i},\bm{q}^{2}_{j},\bm{q}^{2})}\right)\Omega(\bm{q}^{2},m_{\ell})$
$\displaystyle\hskip
62.59596pt\left.+\,2\,\left(m_{\ell}^{2}-\frac{\bm{q}_{i}^{2}+\bm{q}_{j}^{2}-\bm{q}^{2}}{4}-\frac{\bm{q}_{i}^{2}\,\bm{q}_{j}^{2}\,\bm{q}^{2}}{\lambda(\bm{q}^{2}_{i},\bm{q}^{2}_{j},\bm{q}^{2})}\right)\text{C}_{0}\left(q^{2}_{i},q^{2}_{j},q^{2},m_{\ell},m_{\ell},m_{\ell}\right)\right\\}~{},$
$\displaystyle\Omega(x,m_{\ell})$ $\displaystyle\equiv$
$\displaystyle\sqrt{x\,(x+4\,m^{2}_{\ell})}\,\ln\left[\frac{2\,m^{2}_{\ell}+x+\sqrt{x\,(x+4\,m^{2}_{\ell})}}{2\,m^{2}_{\ell}}\right]~{}.$
(23)

The Kallen function is denoted by
$\lambda(\bm{q}^{2}_{i},\bm{q}^{2}_{j},\bm{q}^{2})=\bm{q}^{4}_{i}+\bm{q}^{4}_{j}+\bm{q}^{4}-2\,\bm{q}^{2}_{i}\,\bm{q}^{2}_{j}-2\,\bm{q}^{2}_{i}\,\bm{q}^{2}-2\,\bm{q}^{2}_{j}\,\bm{q}^{2}$,
$\text{C}_{0}(q^{2}_{i},q^{2}_{j},q^{2},m_{\ell},m_{\ell},m_{\ell})$ is the
3-point scalar Passarino-Veltman loop function, and we define $V_{\ell}$ as a
dimensionful constant,

$V_{\ell}\equiv-\frac{8\,\alpha_{\rm em}^{2}}{3\,m_{\ell}}~{}.$ (24)

This is the value of $V_{\alpha}$ in the heavy-lepton limit [50]. The
effective “potential” induced by these two-nucleon interactions is given by
the Fourier transform of the amplitude,

$\tilde{V}_{ij}=-\int\frac{d\bm{q}_{i}}{(2\pi)^{3}}\int\frac{d\bm{q}_{j}}{(2\pi)^{3}}\,e^{-i\bm{q}_{i}\cdot\bm{x}_{i}-i\bm{q}_{j}\cdot\bm{x}_{j}}(2\pi)^{3}\,\delta(\bm{q}-\bm{q}_{i}-\bm{q}_{j})\,{\cal
M}_{ij}(\bm{q}_{i},\bm{q}_{j})~{},$ (25)

The total contribution to the two-body amplitude is given by the sum over all
(identical, i.e., $V_{ij}\equiv V$ for all pairs) proton pairs in the nucleus,
accounting for the nuclear density by means of the nuclear form factor [118,
97],

$\displaystyle{\cal M}_{S}^{(2)}$ $\displaystyle=$
$\displaystyle-\frac{Z(Z-1)}{2}\int d\bm{x}_{i}\int
d\bm{x}_{j}\,\tilde{V}_{ij}\,\rho(\bm{x}_{i},\bm{x}_{j})$ (26)
$\displaystyle=$ $\displaystyle-\frac{Z(Z-1)}{2}\int d\bm{x}_{i}\int
d\bm{x}_{j}\,\tilde{V}_{ij}\,\int\frac{d\tilde{\bm{q}}_{i}}{(2\pi)^{3}}\int\frac{d\tilde{\bm{q}}_{j}}{(2\pi)^{3}}e^{-i\tilde{\bm{q}}_{i}\cdot\bm{x}_{i}-i\tilde{\bm{q}}_{j}\cdot\bm{x}_{j}}\,F^{(2)}(\tilde{\bm{q}}_{i},\tilde{\bm{q}}_{j})$
$\displaystyle=$
$\displaystyle\frac{Z(Z-1)}{2}\,\int\frac{d\bm{k}}{(2\pi)^{3}}\,{\cal
M}_{S,ij}(\bm{q}_{i},\bm{q}_{j})\,F^{(2)}(\bm{q}_{i},\bm{q}_{j})~{},$

where $Z(Z-1)/2$ is the number of proton pairs and
$F^{(2)}(\bm{q}_{i},\bm{q}_{j})$ is the Fourier transform of the two-proton
nuclear density matrix. In what follows we neglect possible (but unknown)
correlations within proton pairs and express $F^{(2)}$ in terms of the one-
proton charge form factor,

$F^{(2)}(\bm{q}_{i},\bm{q}_{j})\simeq
F^{(1)}(\bar{q}_{i})\,F^{(1)}(\bar{q}_{j})=e^{-\bar{q}_{i}^{2}}\,e^{-\bar{q}_{j}^{2}}\hskip
14.22636pt,\hskip 14.22636pt\bar{q}_{i}\equiv|\bm{q}_{i}|/Q_{0}~{}.$ (27)

Figure 3: The scalar form factor, $\bm{{\cal F}_{S}^{\rm
2b}(t;m_{p},m_{\ell})}$, for two-body interactions, as a function of the
square of the exchanged momentum, $\bm{q}^{2}\simeq-t$. We show results for
four different nuclei: helium (left upper panel), oxygen (right upper panel),
iron (left lower panel) and xenon (right lower panel). Note that, although the
factor $Z\,(Z-1)/2$ is not included, the result depends on the target. The
gray shaded area indicates the range of momentum transfer relevant for
scattering processes of local galactic DM off standard matter.

Figure 4: Contributions to the DM-nucleus scalar amplitude, as a function of
the squared of the exchanged momentum $\bm{q}^{2}\simeq-t$: one-body (black
curves), ${\cal M}_{S}^{\rm 1b}$ and two-body (blue curves), ${\cal
M}_{S}^{\rm 2b}$. The total amplitude is also depicted (red curves). We show
results for four different nuclei: helium (left upper panel), oxygen (right
upper panel), iron (left lower panel) and xenon (right lower panel). The gray
shaded area indicates the range of momentum transfer relevant for scattering
processes of local galactic DM off standard matter.

Using Eq. (22), the two-body amplitude can be written as

$\displaystyle{\cal M}_{S}^{(2)}(\bm{q})$ $\displaystyle=$
$\displaystyle-\frac{Z(Z-1)}{2}\,G\,g_{S}\sum_{\ell=e,\mu,\tau}\int\frac{d\bm{k}}{(2\pi)^{3}}\,\frac{\left(\bm{k}-\frac{\bm{q}}{2}\right)\cdot\left(\bm{k}+\frac{\bm{q}}{2}\right)}{\left(\bm{k}-\frac{\bm{q}}{2}\right)^{2}\,\left(\bm{k}+\frac{\bm{q}}{2}\right)^{2}}\,V_{\alpha}(-\bm{k}+\bm{q}/2,\bm{k}+\bm{q}/2,m_{\ell})\,F^{(2)}(-\bm{k}+\bm{q}/2,\bm{k}+\bm{q}/2)$
(28) $\displaystyle\equiv$
$\displaystyle-\frac{Z(Z-1)}{2}\,G\,g_{S}\,\sum_{\ell=e,\mu,\tau}V_{\ell}\
F_{pp}(\bar{q},m_{\ell})\equiv\frac{Z(Z-1)}{2}\,G\,g_{S}\,\sum_{\ell=e,\mu,\tau}{\cal
F}_{S}^{\rm 2b}(\bar{q},m_{\ell})~{},$

where we have made the substitution $\bm{q}_{i}=-\bm{k}+\bm{q}/2$ and
$\bm{q}_{j}=\bm{k}+\bm{q}/2$, $F_{pp}(\bar{q})$ is defined as the two-proton
form factor in the heavy-lepton limit, and ${\cal F}_{S}^{\rm 2b}$ is the two-
body scalar form factor, which is shown in Fig. 3 for different nuclei.

For ease of understanding the $m_{\ell}$ dependence and for comparison with
previous works, we also compute explicitly the $\bm{q}^{2}\to 0$ limit,

$V_{\alpha}(\bm{q}_{i},\bm{q}_{j},m_{\ell})|_{\bm{q}^{2}\to
0}=V_{\ell}\,\frac{6\,m_{\ell}^{2}}{\bm{k}^{2}}\,\left(1-\frac{2\,m_{\ell}^{2}\,\Omega(\bm{k}^{2},m_{\ell})}{\bm{k}^{2}\,(\bm{k}^{2}+4\,m_{\ell}^{2})}\right)~{}.$
(29)

Note that the heavy-lepton limit, $V_{\alpha}\simeq V_{\ell}$, is easily
recovered. Using Eq. (29), one can see that ${\cal F}_{S}^{\rm
2b}(0,m_{\ell})$ is a decreasing function with $m_{\ell}$, with well-defined
limits. In the heavy-lepton limit, ${\cal F}_{S}^{\rm 2b}(0,m_{\ell}\gg
Q_{0})\simeq-\frac{Q_{0}}{4\,\pi\sqrt{2\,\pi}}\,V_{\ell}$, whereas in the
opposite limit, ${\cal F}_{S}^{\rm 2b}(0,m_{\ell}\ll
Q_{0})\simeq-\frac{3\,m_{\ell}}{16}\,V_{\ell}=\alpha_{\rm em}^{2}/2$ tends to
a constant. Although we do not have a close form for the amplitude in the
$\bm{q}^{2}\to 0$, using these limiting values, it is very well approximated
(using two parameters and with an accuracy better than 0.5%) by

$\displaystyle{\cal M}_{S}^{(2)}\left({\bm{q}^{2}\to 0}\right)$
$\displaystyle=$
$\displaystyle\frac{Z(Z-1)}{2}\,G\,g_{S}\,\sum_{\ell=e,\mu,\tau}{\cal
F}_{S}^{\rm 2b}(0,m_{\ell})$ (30) $\displaystyle\simeq$
$\displaystyle\frac{Z(Z-1)}{2}\,G\,g_{S}\,\sum_{\ell=e,\mu,\tau}\frac{3\,m_{\ell}}{16}\,V_{\ell}\,\frac{-1-1.672\,\left(\frac{m_{\ell}}{Q_{0}}\right)}{1+5.589\,\left(\frac{m_{\ell}}{Q_{0}}\right)+1.672\,\frac{3\,\pi^{3/2}}{2\,\sqrt{2}}\,\left(\frac{m_{\ell}}{Q_{0}}\right)^{2}}~{}.$

Thus, the largest contribution at $\bm{q}^{2}=0$ to two-body interactions
comes from the diagram with the electron loop. This is clearly seen in all
panels of Fig. 3, where we depict the two-body form factor for four different
nuclei as a function of $\bm{q}^{2}$. We have numerically evaluated these form
factors with the use of the publicly available packages Package-X [119] and
Collier [120, 121, 122, 123]. As a consequence of the $m_{\ell}/Q_{0}$
dependence, for light nuclei (large $Q_{0}$), the contribution from the
diagram with the muon loop is more important in relative terms than for
heavier nuclei. Similarly to the one-body contribution, the form factor
decreases with $\bm{q}^{2}$ and in the relevant range of momentum transfer,
the diagram with the electron loop gets very suppressed. Indeed, except for
very heavy targets, the diagram with the muon loop is the most important one
in the range $\bm{q}^{2}\sim\left(10^{-4}-10^{-1}\right)~{}{\rm GeV}^{2}$.
Nevertheless, unlike what happens for one-body interactions (c.f., Figs. 2 and
3), the two-body contribution in that range of momentum transfer is
significantly smaller than at $\bm{q}^{2}=0$. This can be understood from two
features. On one hand, at $\bm{q}^{2}=0$, the main contribution to the two-
body amplitude comes from a single diagram, which is subdominant at larger
$\bm{q}^{2}$. On another hand, the drop in the two-body form factor occurs at
smaller values of the exchanged momentum than in the one-body case due to the
term $F^{(2)}(\bm{q}_{i},\bm{q}_{j})$.

In Fig. 4, we compare these two contributions to the total amplitude of the
DM-nuclei elastic scattering cross section. From the results above, it is easy
to see that, at $\bm{q}^{2}=0$, the ratio of the two-body to one-body
amplitudes scales as $(Z-1)/2$.444Note that, in the limit $m_{\ell}\gg m_{N}$,
we recover the scaling as $(Z-1)\,Q_{0}/m_{N}$, discussed in Ref. [97]. At
$\bm{q}^{2}\sim(10^{-4}-10^{-2})~{}{\rm GeV}^{2}$ both contributions are
dominated by the diagrams with the muon loop, and the ratio approximately
scales as $(Z-1)\,Q_{0}/m_{\mu}$. Indeed, for light nuclei like helium, the
dominant contribution comes from one-body interactions at all relevant
$\bm{q}^{2}$. Obviously, for scattering off hydrogen only one-body
interactions take place. For heavier nuclei, however, the two-body
contribution becomes relatively larger and for nuclei heavier than oxygen, it
is the dominant one. Thus, within certain leptophilic scenarios considered
here, one-body interactions are particularly important when computing the DM
capture rate by solar nuclei, but two-body interactions become comparable or
even the largest contribution to event rates in terrestrial direct detection
experiments.

Let us finally comment on the two-body amplitude for the pseudo-scalar case.
This has been briefly discussed in the context of DM scattering induced by
$F_{\mu\nu}\tilde{F}_{\mu\nu}$ [97]. Note that this effective operator is
obtained in the heavy-lepton limit within our scenario. In that work, the
dominant contribution was considered to be the one from the coherent
scattering of one of the two photons on the proton charge and the incoherent
spin-dependent scattering of the second photon on the nucleon magnetic moment.
Nevertheless, the contribution from the scattering of the two photons on the
protons charge is of similar order and should be included too. These
contributions are suppressed with respect to the scalar case by factors of
order $\bm{q}/m_{N}$. Here, we just write down the form of the amplitude,
including these two contributions. In analogy with the scalar case, for
pseudo-scalar couplings the two-body amplitude takes the form

${\cal
M}_{P,ij}=\frac{\Delta_{P}^{\mu\nu}(q_{i},q_{j},m_{\ell})}{q^{2}_{i}\,q^{2}_{j}}\,\left(\bar{u}_{i}\gamma_{\mu}u_{i}\,\bar{u}_{j}\Gamma_{\nu}u_{j}\right)~{},$
(31)

where
$\Gamma_{\nu}=\\{Q_{j}\,\gamma_{\nu},i\,\kappa_{j}\,\sigma_{\nu\mu}\,\frac{q^{\mu}}{2m_{N}}\\}$,
with $Q_{j}$ the electric charge (in units of $e$) and $\kappa_{i}$ the
anomalous magnetic moment (in units of Bohr magneton) of nucleon $j$, and

$\Delta_{P}^{\mu\nu}(q_{i},q_{j},m_{\ell})=G\,g_{P}\,128\,\pi^{2}\,\alpha_{\rm
em}^{2}\,\left(i\,m_{\ell}\,\epsilon^{\mu\nu\rho\sigma}\,q_{i,\rho}\,q_{j,\sigma}\right)\,\text{C}_{0}\left(q^{2}_{i},q^{2}_{j},q^{2},m_{\ell},m_{\ell},m_{\ell}\right)~{}.$
(32)

For $\Gamma_{\nu}=i\,\kappa_{j}\,\sigma_{\nu\mu}\,\frac{q^{\mu}}{2m_{N}}$ we
get an identical result to that in Ref. [97].555Note, however, that $\kappa$
is the nucleon anomalous magnetic moment, not the nucleon magnetic moment as
quoted in Ref. [97]. In the non-relativistic limit, the two-body amplitude can
be written as

${\cal M}_{P}^{(2{\rm
em})}(\bar{q})=G\,g_{P}\,\sum_{i<j}\,\sum_{\ell=e,\mu,\tau}\int\frac{d\bm{k}}{(2\pi)^{3}}\,\frac{1}{\left(\bm{k}-\frac{\bm{q}}{2}\right)^{2}\,\left(\bm{k}+\frac{\bm{q}}{2}\right)^{2}}\,V_{P}^{\rm
em}(\bm{k},\bm{q})\,S^{(2)}(-\bm{k}+\bm{q}/2,\bm{k}+\bm{q}/2)~{},$ (33)

with

$V_{P}^{\rm em}(\bm{k},\bm{q})=128\,\pi^{2}\,\alpha_{\rm
em}^{2}\,m_{\ell}\,\left(\left[\bm{k}^{2}\,\bm{q}-(\bm{k}\cdot\bm{q})\,\bm{k}\right]\cdot\vec{\mu}_{ij}^{\,+}+\left[\bm{q}^{2}\,\bm{k}-(\bm{k}\cdot\bm{q})\,\bm{q}\right]\cdot\frac{\vec{\mu}_{ij}^{\,-}}{2}\right)\,\text{C}_{0}\left(q^{2}_{i},q^{2}_{j},q^{2},m_{\ell},m_{\ell},m_{\ell}\right)~{},$
(34)

where
$\vec{\mu}_{ij}^{\,\pm}=\bm{s}_{i}\,\kappa_{i}\,Q_{j}/m_{i}\pm\bm{s}_{j}\,\kappa_{j}\,Q_{i}/m_{j}$,
with $\bm{s}_{i}=\bm{\sigma}/2$ the spin of nucleon $i$. Analogously to the
scalar case, $S^{(2)}$ accounts for the two-nucleon spin distribution. On
another hand, for $\Gamma_{\nu}=Q_{j}\,\gamma_{\nu}$, we obtain a very similar
expression,

${\cal M}_{P}^{(2{\rm
ee})}(\bar{q})=G\,g_{P}\,\sum_{i<j}\,\sum_{\ell=e,\mu,\tau}\int\frac{d\bm{k}}{(2\pi)^{3}}\,\frac{1}{\left(\bm{k}-\frac{\bm{q}}{2}\right)^{2}\,\left(\bm{k}+\frac{\bm{q}}{2}\right)^{2}}\,V_{P}^{\rm
ee}(\bm{k},\bm{q})\,S^{(2)}(-\bm{k}+\bm{q}/2,\bm{k}+\bm{q}/2)~{},$ (35)

with

$V_{P}^{\rm ee}(\bm{k},\bm{q})=128\,\pi^{2}\,\alpha_{\rm
em}^{2}\,m_{\ell}\,\left(\left[\bm{k}^{2}\,\bm{q}-(\bm{k}\cdot\bm{q})\,\bm{k}\right]\cdot\vec{\varepsilon}_{ij}^{\,+}+\left[\bm{q}^{2}\,\bm{k}-(\bm{k}\cdot\bm{q})\,\bm{q}\right]\cdot\frac{\vec{\varepsilon}_{ij}^{\,-}}{2}\right)\,\text{C}_{0}\left(q^{2}_{i},q^{2}_{j},q^{2},m_{\ell},m_{\ell},m_{\ell}\right)~{},$
(36)

where
$\vec{\varepsilon}_{ij}^{\,\pm}=\left(\bm{s}_{i}/m_{i}\pm\bm{s}_{j}/m_{j}\right)\,Q_{i}\,Q_{j}$.
Note that, as for scalar interactions, only protons contribute to ${\cal
M}_{P}^{(2{\rm ee})}$.

The leading contribution to the two-body scattering amplitude in the pseudo-
scalar case is the sum of the two terms above and can be written as

${\cal M}_{P}^{(2)}(\bar{q})={\cal M}_{P}^{(2{\rm em})}(\bar{q})+{\cal
M}_{P}^{(2{\rm ee})}(\bar{q})=G\,g_{P}\,\sum_{\ell=e,\mu,\tau}{\cal
F}_{P}^{\rm 2b}(\bar{q},m_{\ell})~{}.$ (37)

Notice that the total amplitude has the same form of Eq. (33), but
substituting the anomalous magnetic moment by the magnetic moment. The ratio
of the one-body to two-body contributions is expected to go as in the scalar
case [97], so the amplitude is expected to be dominated by the former for
interactions off light nuclei and by the latter for interactions off heavy
nuclei. We do not explicitly compute numerically this contribution.

##  III Comparison with previous results and discussion

In the literature, there exist various estimates of the two-photon exchange
diagrams in Fig. 1 (with a lepton triangle loop and two exchanged photons), in
particular, within the context of leptophilic DM scenarios [50, 80]. In those
scenarios, the DM particle would only couple to leptons at tree level, but not
to quarks. Therefore, the scattering off nuclei could only take place via
loop-suppressed processes. If the mediator particle couples to charged leptons
via scalar or pseudo-scalar couplings, the loop diagrams here discussed would
be the dominant ones in DM-nuclei scattering, and could determine the capture
rate in astrophysical objects (like the Sun or the Earth) and the event rates
at direct detection experiments. In this work, we have computed the elastic
scattering amplitudes from one- and two-body processes with two-photon
exchange. We have kept all mass dependencies within the loops and, in the one-
body case, we have obtained analytical results in terms of GPLs. Now, we
compare our results to estimates in the literature and to different
approximations.

A first attempt to evaluate the DM-nucleus scattering cross section in the
case of scalar DM-lepton interactions only considered the two-loop diagrams
shown in Fig. 1 [50], but with the entire nucleus as the target. An effective
operator product expansion was used to compute the two-loop diagram, which was
approximated by first integrating out the lepton triangle-loop in the heavy-
lepton limit, then matching it onto the local Rayleigh operator and finally,
attaching the nucleus electromagnetic current. Thus, integrating out the
lepton loop results in the effective (Rayleigh-like) interaction with photons

$\mathcal{L}_{\rm Ray}=\frac{8\,\alpha^{2}_{\rm
em}}{3\,e^{2}}\,\sum_{\ell=e,\mu,\tau}\frac{G\,g_{S}}{4\,m_{\ell}}\,\Gamma_{\chi}F_{\mu\nu}F^{\mu\nu}=-G\,g_{S}\,\sum_{\ell=e,\mu,\tau}V_{\ell}\,\left(\Gamma_{\chi}\,\frac{1}{4}F_{\mu\nu}F^{\mu\nu}\right)~{},$
(38)

where $G\,g_{S}$ is the dimensionful effective coupling of DM particles to
leptons, and $\Gamma_{\chi}$ is the DM bilinear, with
$\Gamma_{\chi}=\\{\chi^{\dagger}\chi\\}$ (scalar) or
$\\{\bar{\chi}\chi,\bar{\chi}\gamma_{5}\chi\\}$ (fermion). We have factored
out the coupling $g_{\beta}$ of the mediator to leptons, so that we can
reabsorb it into $V_{\ell}$.

In the long wavelength limit (small momentum transfer), leptophilic DM would
couple to nuclei through two photons via coherent one- and two-body666Multi-
body interactions with more than two exchanged photons would also contribute,
although at a suppressed level. interactions [97]. Whereas the amplitude of
one-particle interactions is proportional to the number of nucleons, the one
for two-particle interactions is proportional to the number of nucleon pairs,
and hence to $Z(Z-1)/2$. In practice, the latter acts as a coupling to the
total electromagnetic current777At low momentum transfer and at order
$\mathcal{O}(q^{2}/2m_{N})$, this interaction is only sensitive to the total
electric charge of the nucleus, $Ze$, whereas corrections from the nuclear
magnetic moment are suppressed by the nucleus mass, $q/m_{A}$. of the nucleus
[94]. Therefore, attaching the two-photon interaction to the nucleus
electromagnetic current effectively accounts for two-body interactions [97].
In the heavy nucleus limit, Ref. [50] obtained the amplitude

$\mathcal{M}^{(2)}_{\rm
Ray}=-G\,g_{S}\,Z^{2}\,\sum_{\ell=e,\mu,\tau}\frac{V_{\ell}}{4}\,\frac{2\,\pi^{2}|\bm{q}|}{16\,\pi^{2}}\,\tilde{F}(\bar{q})\,\left(\tilde{\Gamma}_{\chi}\bar{u}_{A}u_{A}\right)~{},$
(39)

where $\tilde{\Gamma}_{\chi}=\\{1,\bar{u}_{\chi}u_{\chi}\\}$ for {scalar,
fermion} DM, and $A$ and $Z$ are the mass and atomic numbers, respectively.
The form factor $\tilde{F}(\bar{q})$ was not computed and its $\bm{q}^{2}$
dependence was neglected by setting $\tilde{F}(\bar{q})=1$. The expression
above obviously vanishes at $\bm{q}^{2}=0$. Nevertheless, it should not be so.
Starting from Eq. (38), and using an effective theory approach for the
nucleus, the amplitude of the one-loop contribution with the two photons
attached to the entire nucleus reads [94, 93, 97]

$\displaystyle\mathcal{M}^{(2)}_{\rm Ray}(\bar{q})$ $\displaystyle=$
$\displaystyle-G\,g_{S}
\,Z^{2}\,\sum_{\ell=e,\mu,\tau}V_{\ell}\,\frac{Q_{0}}{4\,\pi\,\sqrt{2\,\pi}}\,F_{\rm
Ray}(\bar{q})\,\left(\tilde{\Gamma}_{\chi}\bar{u}_{A}u_{A}\right)\equiv-G\,g_{S}\,Z^{2}\,\sum_{\ell=e,\mu,\tau}V_{\ell}\,F^{\rm
eff}_{pp}(\bar{q})\,\left(\tilde{\Gamma}_{\chi}\bar{u}_{A}u_{A}\right)~{},$
(40) $\displaystyle F_{\rm Ray}(\bar{q})$ $\displaystyle=$
$\displaystyle-\sqrt{\frac{2}{\pi}}\,\int^{1}_{0}dx\int^{\infty}_{0}dl\,\frac{l^{2}}{\left(l^{2}+\bar{q}^{2}\,(1-x)\,x\right)^{2}}\,\exp\left[-2\,\left(l^{2}-\bar{q}^{2}\,\left(x\,(1-x)-\tfrac{1}{2}\right)\right)\right]$
(41)
$\displaystyle\times\Big{\\{}\cosh\left[2\,l\,\bar{q}\,(1-2\,x)\right]-\frac{l^{2}-\bar{q}^{2}\,(1-x)\,x+\tfrac{1}{2}}{l\,\bar{q}\,(1-2\,x)}\sinh\left[2\,l\,\bar{q}\,\left(1-2\,x\right)\right]\Big{\\}}~{},$

where $F_{\rm Ray}$ is defined so that $F_{\rm Ray}(0)=1$. Up to order ${\cal
O}(\bar{q}^{3})$, the two-proton form factor in this limit is given by [97]

$F_{pp}^{\rm
eff}(\bar{q})\simeq\frac{Q_{0}}{4\,\pi\sqrt{2\,\pi}}\,e^{-\bar{q}^{2}/2}\,\left[1-\frac{\pi^{3/2}}{2\,\sqrt{2}}\,\bar{q}+\frac{5}{3}\,\bar{q}^{2}\right]~{}.$
(42)

This is the result for $F_{pp}$ we also obtain in the limit when the Wilson
coefficient corresponding to the triangle loop is taken to be momentum-
independent [94, 97], as is the case of the heavy-lepton limit,
$m_{\ell}\gtrsim Q_{0}$, for which $V_{\alpha}=V_{\ell}$ [50]. The ratio of
our result to the above estimate is given by

$R^{(2)}(\bar{q})=\frac{{\cal M}_{S}^{(2)}(\bar{q})}{{\cal M}^{(2)}_{\rm
Ray}(\bar{q})}=-\frac{Z-1}{2\,Z}\,\frac{\sum_{\ell=e,\mu,\tau}{\cal
F}_{S}^{\rm 2b}(\bar{q},m_{\ell})}{\sum_{\ell=e,\mu,\tau}V_{\ell}\,F^{\rm
eff}_{pp}(\bar{q})}~{}.$ (43)

Let us first focus on the $\bm{q}^{2}=0$ case. For $m_{\ell}\gtrsim Q_{0}$,
the heavy-lepton limit is a very good approximation for tau and even for muon
leptons in the loop. However, it significantly overestimates the (absolute)
value of ${\cal M}^{(2)}$ with electrons running in the loop. Note that this
is particularly relevant, as the contribution from the loop with the lightest
lepton is the dominant one. For $m_{\ell}\gtrsim Q_{0}$, our two-body
amplitude, Eq. (30), is smaller than Eq. (40) by a factor $(Z-1)/(2\,Z)$.
Nevertheless, for $m_{\ell}\ll Q_{0}$,

$R^{(2)}_{\ell}(0)=\frac{{\cal M}_{S,\ell}^{(2)}(0)}{{\cal M}^{(2)}_{{\rm
Ray},\ell}(0)}\,\stackrel{{\scriptstyle{m_{\ell}\ll
Q_{0}}}}{{\simeq}}\,\frac{Z-1}{2\,Z}\,\frac{3\,\pi\sqrt{2\,\pi}}{4}\,\frac{m_{\ell}}{Q_{0}}\,\stackrel{{\scriptstyle
m_{\ell}=m_{e}}}{{\sim}}\,(0.003-0.015)~{}.$ (44)

Therefore, given the typical size of nuclei, $Q_{0}\simeq(0.1-0.3)$ GeV, the
heavy-lepton limit overestimates the amplitude by two-three orders of
magnitude (the difference being larger for lighter nuclei) and thus, the cross
section by four-five orders of magnitude. As discussed above, at higher
momenta, the contribution from the diagram with electrons in the loop gets
suppressed and the diagram with muons in the loop becomes the most important
one (see Fig. 2). However, with the approximation in Eq. (40), the relative
importance of the different lepton loops does not change with the value of the
momentum transfer and the diagram with the electron loop is always the
dominant one. As a result, the ratio $R_{\ell}^{(2)}$ decreases with
$\bm{q}^{2}$, and at $\bm{q}^{2}\sim(10^{-3}-10^{-2})~{}{\rm GeV}^{2}$ it is
up to an order of magnitude smaller than at $\bm{q}^{2}=0$. In summary, the
heavy-lepton approximation overestimates the two-body contribution to the
scattering cross section by more than five orders of magnitude in the relevant
range of momentum transfer.

The one-body interaction was not considered for leptophilic DM models in Ref.
[50], but it was the only contribution considered in Ref. [80]. Moreover, it
was accounted for within more general approaches [93, 97], which start from an
effective description of DM interactions with two photons via a Rayleigh-like
operator, $\Gamma_{\chi}F_{\mu\nu}F^{\mu\nu}$. The one-body amplitude was
taken to be $(Z\,f_{p}+(A-Z)\,f_{n})\,F^{(1)}(t/Q_{0}^{2})$, with $f_{p}$ and
$f_{n}$ the proton and neutron couplings. Whereas Ref. [97] followed a purely
phenomenological approach without specifying $f_{p}$ or $f_{n}$, Ref. [93]
considered the contribution to the quark $\Gamma_{\chi}\bar{q}q$ and gluon
$\Gamma_{\chi}G_{\mu\nu}G^{\mu\nu}$ operators, via mixing with the Rayleigh
operator. Then, the matrix element with nucleons was connected to the
amplitudes at the quark and gluon level with the usual procedure at vanishing
momentum transfer. In this paper, we have followed a different approach. Given
that in the scenario under consideration, interactions with nucleons occur via
photon exchange, there is no need to go through couplings to quarks. We have
computed one-body interactions directly using the electromagnetic coupling to
nucleons. In this way, we save one step and avoid having to deal with dressed
quark propagators in the two-loop diagrams in Fig. 1. Our result is given in
Eq. (11), and inserting the external DM particles implies multiplying by $G$.
We compare our result for the scalar mediator case to [93, 80]

${\cal M}^{(1)}_{qq}=(G\,g_{S})\,\left[\sum_{N=p,n}\,{\cal
N}_{N}\,\left(\sum_{q=u,d,s}\frac{m_{N}}{m_{q}}\,f_{q}^{N}\,{\cal
C}_{q}+\frac{2}{27}\,f_{G}^{N}\sum_{Q=c,b,t}\frac{m_{N}}{m_{Q}}\,\Delta_{Q}\,{\cal
C}_{Q}\right)\right]\,F^{(1)}(\bar{q})$ (45)

where ${\cal N}_{p}=Z$ and ${\cal N}_{n}=A-Z$,
$\Delta_{Q}=1+11\,\alpha_{s}(m_{Q})/(4\pi)$ with $\alpha_{s}(\mu)$ the strong
coupling at scale $\mu$, and the assumed values of the hadronic matrix
elements at $t=0$ are
$\\{f_{u}^{p},\,f_{d}^{p},\,f_{s}^{p},\,f_{G}^{p}\\}=\\{0.016,\,0.029,\,0.043,\,0.912\\}$
and
$\\{f_{u}^{n},\,f_{d}^{n},\,f_{s}^{n},\,f_{G}^{n}\\}=\\{0.014,\,0.034,\,0.043,\,0.909\\}$.
For the quark masses we use the values quoted at the Particle Data Group
[124]. The terms $(G\,g_{S})\,{\cal C}_{q}$ are the Wilson coefficients of the
operator $\Gamma_{\chi}\bar{q}q$.

In Ref. [80], the integration of the two-loop diagram seems to be performed in
the limit of vanishing square momentum transfer ($t=0$). Implicitly, it is
claimed that the scalar and pseudo-scalar form factors are equal in that
limit.888In Ref. [80] only fermionic leptophilic DM is considered. In the case
of scalar or pseudo-scalar mediators, DM bilinears including $\gamma_{5}$ are
claimed not to induce couplings to quarks at one or two loops. We do not agree
with this statement. Our results also apply to those cases. Although details
about the followed steps are not provided, we can compare our results with the
approximate analytical expression that is obtained (Eq. (4.14) in Ref. [80]),
which we reproduce here in slightly different notation,

${\cal C}_{q}^{\rm
BBR}\left(z_{q\ell}\equiv\frac{m_{q}}{m_{\ell}}\right)=Q_{q}^{2}\,\left(\frac{\alpha_{\rm
em}^{2}}{2\,\pi^{2}}\right)\,\sum_{\ell=e,\mu,\tau}\,\frac{z_{q\ell}}{1-\left(\frac{z_{q\ell}}{2}\right)^{2}}\,\ln\left(\frac{2}{z_{q\ell}}\right)\hskip
28.45274pt~{},$ (46)

where $Q_{q}$ is the electric charge of quark $q$. On another hand, the Wilson
coefficients obtained from Ref. [93], after computing the effective coupling
of the Rayleigh operator in the heavy-lepton limit ($V_{\alpha}=V_{\ell}$)
within our leptophilic scenario, can be written as

${\cal C}^{\rm EFT}_{q}(z_{q\ell})=-\left(\frac{1}{4\,\pi\,\alpha_{\rm
em}}\,\sum_{\ell=e,\mu,\tau}V_{\ell}\,m_{q}\right)\,3\,Q_{q}^{2}\,\left(\frac{\alpha_{\rm
em}}{\pi}\right)\,\ln\left(\frac{\Lambda^{2}}{m_{y}^{2}}\right)=Q_{q}^{2}\,\left(\frac{\alpha_{\rm
em}^{2}}{2\,\pi^{2}}\right)\,\sum_{\ell=e,\mu,\tau}4\,z_{q\ell}\,\ln\left(\frac{\Lambda^{2}}{m_{y}^{2}}\right)~{},$
(47)

where $m_{y}=m_{N}$ for $q=\\{u,d,s\\}$ and $m_{y}=m_{Q}$ for $Q=\\{c,b,t\\}$,
and we will take $\Lambda=1$ TeV, although the value of the UV scale does not
affect our conclusions. Note that, in addition to neglecting the momentum
dependence of the hadronic matrix elements, these effective couplings are not
valid for $\bm{q}^{2}\gtrsim m_{\ell}^{2}$, so we can only strictly compare
this estimate to our result at very low values of the momentum transfer.
Therefore, we will only compare the above two estimates to our results at
$t=0$. For the scalar case, the ratio of our calculation to the ones above is

$R^{(1)}(0)\equiv\frac{{\cal M}^{(1)}(0)}{{\cal
M}^{(1)}_{qq}(0)}=\frac{Z\,\sum_{\ell=e,\mu,\tau}{\cal F}_{S}^{\rm
1b}(t=0;z_{p\ell})}{\sum_{N=p,n}{\cal
N}_{N}\,\left(\sum_{q=u,d,s}\frac{m_{N}}{m_{q}}\,f_{q}^{N}\,{\cal
C}_{q}+\frac{2}{27}\,f_{G}^{N}\sum_{Q=c,b,t}\frac{m_{N}}{m_{Q}}\,\Delta_{Q}\,{\cal
C}_{Q}\right)}\simeq\begin{cases}1.5&\text{with ${\cal C}_{q}^{\rm
BBR}$}~{},\\\\[8.61108pt] 6\times 10^{-4}&\text{with ${\cal C}^{\rm
EFT}_{q}$}~{}.\\\ \end{cases}$ (48)

Note that if we compute the one-body two-loop interaction at the quark level
and then we connect it with the scattering off nuclei as done above (i.e., we
set ${\cal C}_{q}=\sum_{\ell=e,\mu,\tau}\,{\cal F}_{S}^{\rm
1b}(t=0;m_{q},m_{\ell})$), the ratio we obtain is $R^{(1)}(0)\simeq 0.35$.
This implies that the treatment at the quark level for these type of
interactions overestimates the value of the amplitude (at $t=0$) by a factor
of a few, unless the up and down quark masses were actually larger by a
similar factor.

Clearly, the approach using ${\cal C}^{\rm EFT}_{q}$ overestimates the one-
body contribution to the amplitude by about three orders of magnitude. Like in
the two-body processes, the discrepancy has to do with the heavy-lepton limit
being not appropriate for the dominant diagram (the one with the electron
loop) and largely overestimating the correct result.

Let us now briefly discuss the similar results obtained with our approach at
the nucleon level and the estimate using ${\cal C}_{q}^{\rm BBR}$ at the quark
level. First, note that when comparing the results using ${\cal C}^{\rm BBR}$
at the nucleon and quark levels, a similar ratio ($\sim 0.2$) as the one
within our approach is obtained. Therefore, comparing the result in Ref. [80]
to ours, both either at quark or nucleon level, the ratio is $\sim 4.4$ and
$\sim 7.6$, respectively. Nevertheless, when comparing our result at nucleon
level to that of Ref. [80] at quark level, numerically the sum of all one-body
diagrams is similar to our result, Eq. (48). Their relative contributions,
however, are different. Whereas in our case the contributions with the
electron and muon loop are similar and that with the tau loop is about a
factor of three smaller, using Eq. (46) the electron loop is much larger than
the other two contributions. The relative importance of diagrams with a
different lepton in the loop can be understood from the calculation at the
quark or nucleon level (i.e., depending on the internal propagator). Note that
when comparing the calculations at nucleon and at quark level, one should
compare ${\cal F}_{S}^{\rm 1b}(z_{p\ell})$ to $(z_{p\ell}/z_{q\ell})\,{\cal
C}_{q}(z_{q\ell})$ (or to ${\cal F}_{S}^{\rm 1b}(z_{q\ell})\,(z_{q\ell})$).
The former is a growing function with $z_{p\ell}$ which flattens at
$z_{p\ell}\gtrsim 1$ and the latter is a decreasing function with $z_{q\ell}$.
This explains the relative importance of the three diagrams in each case.

Nevertheless, the major difference between our one-body form factor and that
in Ref. [80] is the result of the two-loop integral itself, Eqs. (17) and
(46). For lepton masses in the loop larger than the mass of the external
fermion, $z\ll 1$, Eq. (46) is smaller than Eq. (17) by a factor of $\sim 2$,

$\frac{{\cal F}_{S}^{\rm 1b}(t=0;z\ll 1)}{{\cal C}^{\rm BBR}(z\ll 1)}\simeq
2\,\left(\frac{\frac{13}{12}-\ln(z)}{\ln 2-\ln(z)}\right)~{}.$ (49)

This limit applies to the diagrams with the muon and tau loop, for light
external quarks, which are the ones that contribute the most. In the opposite
limit, $z\gg 1$, the disagreement is larger and grows with $z$,

$\frac{{\cal F}_{S}^{\rm 1b}(t=0;z\gg 1)}{{\cal C}^{\rm BBR}(z\gg
1)}\simeq\frac{\pi^{2}}{8}\,\frac{z}{\ln\left(z/2\right)}~{}.$ (50)

This limit always applies to the diagram with the electron loop, but depending
on the mass of the external fermion, it might or not apply to the other two
diagrams. In particular, for light external quarks, it does not apply to the
diagrams with the muon or tau loop. On the contrary, it is a good
approximation for all diagrams when the external fermion is considered to be a
nucleon. In summary, it is evident that these results are quite different,
mainly for $z\gtrsim 1$ and it is a bit surprising the contributions to the
amplitude in our case and in Ref. [80], when summed over all diagrams, come up
being similar.

Finally, let us mention that in Ref. [50] the pseudo-scalar form factor is
claimed to be zero at all orders. This is also a consequence of integrating
out first the lepton loop. By doing so, other unitarity cuts contributing to
the amplitude are neglected [31, 33, 125]. As we have shown by performing the
one-body two-loop (and have outlined for the two-body one-loop) calculation,
the pseudo-scalar form factor is not zero, and can actually be larger than the
scalar one. We do not further comment on the contribution from two-body
interactions in this case, which is proportional to $Z$ and to an spin-
dependent part (see also Ref. [97] for a discussion).

##  IV Conclusions

In this work, we have considered the scattering of dark matter (DM) off
standard matter nuclei, within leptophilic scenarios, by considering loop-
induced interactions via a scalar and a pseudo-scalar mediator. For these
cases, the leading coupling to nucleons would arise from Feynman diagrams that
contain a lepton triangle loop and two photons exchanged with the external
nucleons, either via one-body two-loop or via two-body one-loop interactions.
For some kinematic regimes, the considered loop-induced processes could become
the dominant contributions (with respect to tree-level DM-electron
interactions) to the rates in direct detection experiments and to the capture
of DM particles in celestial objects.

We evaluated analytically the form factors coming from one- and two-loop
diagrams containing an internal triangle loop of a (massive) lepton and
external (massive) fermion currents, with different masses, and for arbitrary
value of the momentum transfer. Our analytical calculations show that results
obtained by not including the proper dependence on the kinematic variables may
lead to incorrect conclusions.

In particular, the two-loop form factors have been evaluated analytically by
making use of advanced multi-loop calculus techniques, based on IBPs, and
Magnus exponential method for differential equations. The calculations have
been carried out within the dimensional regularization scheme, and the final
expressions of the form factors have been obtained by taking the smooth four
dimensional limit, verifying the finiteness of the results. We have expressed
the form factors in terms of GPLs and transcendental constants (up to weight
four), by keeping the complete dependence on the lepton mass, on the nucleon
mass, and on the momentum transfer. The two-loop scalar form factor had
already been considered in recent studies involving the Higgs boson decay into
massive quark pairs, with a quark triangle loop and two-gluon exchange [16,
18]. We recomputed it here for leptophilic DM scenarios. The results of the
pseudo-scalar form factors are presented here for the first time. Within
leptophilic DM scenarios, previous estimates of the amplitude of the same two-
loop diagrams had been obtained under different approximations. Two-body one-
loop processes had also been considered, under different approximations, in
effective DM scenarios, where DM couples to photons through a Rayleigh-like
operator [94, 93, 97]. We have evaluated the two-body contributions in the
non-relativistic limit, but keeping the dependence on all the relevant scales
and computing the lepton triangle loop without any further approximation.

Notably, we show that loop-induced one-body and two-body interactions could
both significantly contribute to the scattering of DM particles off nuclei in
certain leptophilic scenarios. For one-body interactions, the largest
contributions at small momentum transfer come from the diagrams with the
electron and muon loops. For two-body processes, the dominant diagram at small
momenta is that with the electron loop. At larger momenta, in particular at
momenta relevant for local galactic DM interactions,
$\bm{q}^{2}\sim(10^{-4}-10^{-1})~{}{\rm GeV}^{2}$, the diagram with the
electron loop is exponentially suppressed and the dominant contribution is
from the diagram with muon in the loop. The relative importance of the three
different diagrams (with different lepton loops) varies as a function of the
momentum transfer, but also depends on several other factors. Whereas in the
case of one-body interactions the three diagrams are comparable at $t=0$, for
two-body interactions the diagram with the electron loop is significantly
larger than the other two at $t=0$, and is larger for heavier targets. At
higher momentum transfer, the scalar one-body form factor is just slightly
smaller than at $t=0$, so the $t=0$ approximation provides a reasonable
result. For two-body interactions, the diagram with the tau loop is always
subdominant at all relevant momentum transfer. The differences in relative
importance between the other two diagrams is larger than for one-body
interactions, so using the $t=0$ form factor at other $t$ could result in
overestimating this contribution to the amplitude by about an order of
magnitude.

In previous studies of these type of leptophilic DM scenarios [50, 80], only
one of these contributions was considered and several approximations were
included. We have compared our results to those estimates and conclude that
the predicted estimate for the amplitude of one-body interactions at $t=0$ in
Ref. [80] (performed at the quark level) has approximately the same numerical
value as that of our full two-loop calculation (performed at nucleon level),
although the analytical result is quite different. Thus the previously
estimated scattering cross section for light nuclei, which is dominated by
one-body interactions, is numerically off just by a factor of a few units.
Note, however, that by following the effective approach of Ref. [93], the one-
body amplitude is overestimated by about three orders of magnitude. On another
hand, for nuclei heavier than oxygen, for which the amplitude is dominated by
the two-body contribution, the heavy-lepton approximation [50, 94, 93, 97]
overestimates the scattering cross section by no less than four orders of
magnitude, in the relevant range of momentum transfer. Moreover, we have
showed that the two-loop pseudo-scalar form factor has a non-vanishing
expression, although it had been claimed to be vanishing at all orders [50].

In summary, we have computed and critically examined two-photon exchanged
diagrams, required for the calculation of the DM-nucleus scattering cross
section in some leptophilic DM scenarios. We have presented new analytical
calculations for one-body two-loop and two-body one-loop DM-nucleus
interactions, which are needed when evaluating event rates in direct detection
experiments and in the computation of the DM capture rate by celestial
objects. The calculations in this work can be further extended to consider the
form factors for axial-vector mediators, where anomalous terms might also play
a role. We leave this application to future studies.

## Acknowledgments

We thank Seva Chestnov and Paride Paradisi for clarifying discussions. RG is
supported by MIUR grant PRIN 2017FMJFMW. The work of PM is partially supported
by the Supporting TAlent in ReSearch at Padova University (UniPD STARS Grant
2017 “Diagrammalgebra"). SPR is supported by the Spanish FEDER/MCIU-AEI grant
FPA2017-84543-P, and partially, by the Portuguese FCT (UID/FIS/00777/2019 and
CERN/FIS-PAR/0004/2019). SPR also acknowledges support from the European ITN
project HIDDeN (H2020-MSCA-ITN-2019//860881-HIDDeN). AP is supported by the
Swiss National Science Foundation under grant number 200020-175595. FG is
supported by Fondazione Cassa di Risparmio di Padova e Rovigo (CARIPARO). SPR
and RG also thank the Galileo Galilei Institute for hospitality. Finally, PM
and SPR acknowledge Silvia Pascoli’s encouragement and stimulating discussions
over the years.

##  Appendix A Differential equations for one-body interactions

In this appendix, we provide the details about the calculation of the Master
Integrals (MIs) required for the determination of the one-body form factors
$\mathcal{F}_{S}^{1b}(t;m_{N},m_{\ell})$ and
$\mathcal{F}_{P}^{1b}(t;m_{N},m_{\ell})$. We consider the generic case, with
complete dependence on the variables $\\{t,m_{N},m_{\ell}\\}$, and discuss, as
well, two limiting behaviors: the soft limit (i.e., $t\to 0$), and the equal-
mass limit (i.e., $m_{N}\to m_{\ell}$). In particular, in the soft limit,
$\mathcal{F}_{S}^{1b}(t\to 0;m_{N},m_{\ell})$ is evaluated by setting $t=0$ in
the integrand, which greatly simplifies the integral evaluation; whereas
$\mathcal{F}_{P}^{1b}(t\to 0,m_{N},m_{\ell})$ admits an asymptotic expansion,
performed after integration. The equal-mass limit is studied by verifying that
the expressions derived by expanding $\mathcal{F}_{S}^{1b}(t;m_{N},m_{\ell})$
and $\mathcal{F}_{P}^{1b}(t;m_{N},m_{\ell})$ agree with the ones obtained from
a calculation where the choice $m_{N}=m_{\ell}$ is set ${\it ab\ initio}$ in
the integrands. The evaluation of the MIs for the equal-mass case is
considered separately and the agreement between the two ways of computing the
form factors in this limit provides a check on the expressions obtained for
arbitrary values of the kinematic variables.

###  A.1 Master integrals for the general kinematic case

In this subsection, we summarize the calculation of the MIs required for the
determination of the one-body form factors
$\mathcal{F}_{S}^{1b}(t;m_{N},m_{\ell})$ and
$\mathcal{F}_{P}^{1b}(t;m_{N},m_{\ell})$ through the differential equations
method along the lines of Ref. [16]. As discussed in the main text, all the
integrals of the type of Eq. (14) can be reduced via IBPs to linear
combinations of 20 independent MIs, which obey coupled first order
differential equations in the kinematic variables $x$ and $y$ introduced in
Eq. (15). We do this using the IBPs code Reduze2 [100] and LiteRed [110]. The
two systems of differential equations of Eq. (II.1) can be combined into the
total differential of the vector $\mathbf{I}$ of basis MIs,

$d\,\mathbf{I}(\epsilon,x,y)=\epsilon\,d\mathbb{A}(x,y)\,\mathbf{I}(\epsilon,x,y)\qquad\text{with}\qquad
df=dx\frac{\partial f}{\partial x}+dy\frac{\partial f}{\partial y}\ .$ (51)

The factorization of the dimensional regulator $\epsilon$ from the kinematics
dependence of Eq. (51), which is encoded in the matrix $\mathbb{A}_{j}$,
greatly simplifies the task of determining the expression of the MIs as a
Laurent series around $\epsilon=0$ (i.e., $d=4$), since it decouples the
solution of the differential equations order by order in $\epsilon$. This
factorization is, however, not a general property of the differential
equations for Feynman integrals (which usually exhibit a rational dependence
on $\epsilon$) and it is rather a consequence of a specific choice of the MIs.
In fact, for a given family of Feynman integrals, the basis of MIs is not
unique, as any complete set of integrals, which are linearly independent under
IBPs, can be taken as a valid basis. Although different basis of MIs are
completely equivalent, only certain ones fulfill canonical
$\epsilon$-factorized differential equations of the type of Eq. (51). For the
problem under consideration, the following basis has been identified through
the Magnus-Dyson algorithm [111],

$\displaystyle\text{I}_{1}$
$\displaystyle=\epsilon^{2}\,\mathcal{T}_{1}~{},\qquad$
$\displaystyle\text{I}_{2}=\epsilon^{2}\,\mathcal{T}_{2}~{},$
$\displaystyle\text{I}_{3}$
$\displaystyle=\epsilon^{2}\,\lambda_{\ell}\,\mathcal{T}_{3}~{},\qquad$
$\displaystyle\text{I}_{4}=-\epsilon^{2}\,t\,\mathcal{T}_{4}~{},$
$\displaystyle\text{I}_{5}$
$\displaystyle=\epsilon^{2}\,\left(\frac{1}{2}\,\left(-t+\lambda_{\ell}\right)\,\mathcal{T}_{4}+\lambda_{\ell}\,\mathcal{T}_{5}\right)~{},\qquad$
$\displaystyle\text{I}_{6}=-\epsilon^{2}\,t\,\mathcal{T}_{6}~{},$
$\displaystyle\text{I}_{7}$
$\displaystyle=\frac{\epsilon^{2}\,m_{N}^{2}\,\left(t+\lambda_{\ell}\right)\,\rho_{N}}{\left(\lambda_{N}+t\right)\,\rho_{\ell}}\,\left(\mathcal{T}_{7}+2\,\mathcal{T}_{8}\right)~{},\qquad$
$\displaystyle\text{I}_{8}=\epsilon^{2}\,m_{N}^{2}\,\mathcal{T}_{8}~{},$
$\displaystyle\text{I}_{9}$
$\displaystyle=\epsilon^{2}\,\lambda_{\ell}\,\mathcal{T}_{9}~{},\qquad$
$\displaystyle\text{I}_{10}=-\epsilon^{2}\,t\,\lambda_{\ell}\,\mathcal{T}_{10}~{},$
$\displaystyle\text{I}_{11}$
$\displaystyle=\epsilon^{3}\,\lambda_{N}\,\mathcal{T}_{11}~{},\qquad$
$\displaystyle\text{I}_{12}=\frac{\epsilon^{2}\,\lambda_{\ell}}{4\,t}\,\Big{(}\left(t-\lambda_{N}\right)\,\left(\mathcal{T}_{4}+2\,\mathcal{T}_{5}\right)-4\,m_{N}^{2}\,\lambda_{N}\,\mathcal{T}_{12}\Big{)}~{},$
$\displaystyle\text{I}_{13}$
$\displaystyle=\epsilon^{3}\,\lambda_{N}\,\mathcal{T}_{13}~{},\qquad$
$\displaystyle\text{I}_{14}=\epsilon^{3}\,(-1+2\,\epsilon)\,t\,\mathcal{T}_{14}~{},$
$\displaystyle\text{I}_{15}$
$\displaystyle=\epsilon^{3}\,\lambda_{\ell}\,\lambda_{N}\,\mathcal{T}_{15}~{},\qquad$
$\displaystyle\text{I}_{16}=\epsilon^{3}\,\lambda_{N}\,\mathcal{T}_{16}~{},$
$\displaystyle\text{I}_{17}$
$\displaystyle=\epsilon^{3}\,\lambda_{N}\,\mathcal{T}_{17}~{},\qquad$
$\displaystyle\text{I}_{18}=\frac{\epsilon^{2}}{t}\,\Big{(}\lambda_{\ell}\,\left(-t+\lambda_{N}\right)\,\mathcal{T}_{9}+\epsilon\,\left(t-\lambda_{\ell}\right)\,\lambda_{N}\,\mathcal{T}_{17}+(-1+2\,\epsilon)\,\lambda_{\ell}\,\lambda_{N}\,\mathcal{T}_{18}\Big{)}~{},$
$\displaystyle\text{I}_{19}$
$\displaystyle=\frac{\epsilon^{2}}{2\,t}\,\Big{(}t\,\left(\lambda_{\ell}-t\right)\,\mathcal{T}_{3}-2\,t\,m_{\ell}^{2}\,\left(\mathcal{T}_{7}+2\,\mathcal{T}_{8}\right)$
$\displaystyle+\left(4\,t\,m_{\ell}^{2}+\lambda_{\ell}\,\left(\lambda_{N}-t\right)\right)\,\mathcal{T}_{9}+\epsilon\,\frac{4\,t^{2}\,m_{N}^{2}}{\lambda_{N}+t}\,\mathcal{T}_{16}$
$\displaystyle\qquad+\epsilon\,\left(\lambda_{N}\,\left(t-\lambda_{\ell}\right)-4\,t\,m_{\ell}^{2}\right)\,\,\mathcal{T}_{17}\,\,\,$
$\displaystyle+(2\,\epsilon-1)\,\left(4\,t\,m_{\ell}^{2}+\lambda_{\ell}\,\lambda_{N}-t^{2}\right)\,\mathcal{T}_{18}+2\,t^{2}\,\left(m_{\ell}^{2}-m_{N}^{2}\right)\,\mathcal{T}_{19}\Big{)}~{},$
$\displaystyle\text{I}_{20}$
$\displaystyle=-\epsilon^{4}\,t\,\lambda_{N}\,\mathcal{T}_{20}~{},$ (52)

with $\mathcal{T}_{i}$’s depicted in Fig. 5. In Eqs. (52) we have introduced
the notation

$\lambda_{i}=\sqrt{-t}\,\sqrt{4\,m^{2}_{i}-t}~{},\quad\rho_{i}=\sqrt{\frac{2\,m^{2}_{i}-t-\lambda_{i}}{m^{2}_{i}}},\quad\text{with:
}i=N,\ell~{}.$ (53)

(1)  $\mathcal{T}_{1}$ (2)  $\mathcal{T}_{2}$ (3)  $\mathcal{T}_{3}$ (4)
$\mathcal{T}_{4}$

(5)  $\mathcal{T}_{5}$ (6)  $\mathcal{T}_{6}$

(7)  $\mathcal{T}_{7}$ (8)  $\mathcal{T}_{8}$ (9)  $\mathcal{T}_{9}$ (10)
$\mathcal{T}_{10}$

(11)  $\mathcal{T}_{11}$ (12)  $\mathcal{T}_{12}$

(13)  $\mathcal{T}_{13}$ (14)  $\mathcal{T}_{14}$ (15)  $\mathcal{T}_{15}$
(16)  $\mathcal{T}_{16}$

(17)  $\mathcal{T}_{17}$ (18)  $\mathcal{T}_{18}$

(19)  $\mathcal{T}_{19}$ (20)  $\mathcal{T}_{20}$

Figure 5: Master integrals required for the evaluation of the scalar and
pseudo-scalar form factors, corresponding to the two-loop diagrams in Fig. 1.
In each diagram, we indicate the off-shell external mediator (dashed lines),
fermions with mass $m_{\ell}$ (blue lines) and $m_{N}$ (red lines), massless
photons (black lines) and squared denominators (dots).

A few comments are in order. First, the $\epsilon$ dependence of the
coefficients in Eq. (52) has been adjusted in such a way that all integrals
are finite in the four dimensional limit. Therefore, they admit, rather than a
Laurent series around $\epsilon=0$, a Taylor expansion of the type

$\mathbf{I}(\epsilon,x,y)=\sum_{k=0}^{\infty}\,\mathbf{I}^{(k)}(x,y)\,\epsilon^{k}~{}.$
(54)

Second, we redefine at this stage of the computation the loop integration
measure as

$\int\widetilde{d^{d}k_{i}}=\left(\frac{m^{2}_{l}}{\mu^{2}}\right)^{\epsilon}\,\int\frac{d^{d}k_{i}}{i\,\pi^{d/2}\,\Gamma_{\epsilon}}~{},$
(55)

with $\mu$ being the ’t Hooft scale of dimensional regularization and
$\Gamma_{\epsilon}=\Gamma(1+\epsilon)$ is the $\Gamma$ function. In this way,
all $\text{I}_{i}$ are dimensionless functions of the variables $x$ and $y$
and the two-loop tadpole integral is conveniently normalized to one, namely
$\text{I}_{1}(\epsilon,x,y)=1$. The standard normalization $1/(2\pi)^{d}$ can
be then be restored in the final expression of $\mathcal{F}_{S}^{1b}$ and
$\mathcal{F}_{P}^{1b}$. Finally, it is worth noticing that, while all the $20$
listed MIs appear in the decomposition of the scalar form factor
$\mathcal{F}_{S}^{1b}$, the pseudo-scalar one $\mathcal{F}_{P}^{1b}$ turns out
to be decomposed as a linear combination of only $10$ MIs, namely
$\\{\text{I}_{4},\text{I}_{5},\text{I}_{7},\text{I}_{8},\text{I}_{11},\text{I}_{12},\text{I}_{14},\text{I}_{16},\text{I}_{17},\text{I}_{20}\\}$.

Owing to the $\epsilon$-factorization of the differential equation, it is
straightforward to write a formal expression of the Taylor coefficients
$\mathbf{I}^{(k)}(x,y)$ introduced in Eq. (54) in terms of iterated integrals
over the kernel matrix $\mathbb{A}(x,y)$,

$\mathbf{I}^{(k)}(x,y)=\sum_{i=0}^{k}\,\int_{\gamma}\underbrace{d\mathbb{A}\dots
d\mathbb{A}}_{\text{{\hbox{i}} times}}\,\mathbf{I}^{(k-i)}(x_{0},y_{0})$ (56)

where $\gamma$ is a regular path in the $xy$-plane and
$\mathbf{I}^{(i)}(x_{0},y_{0})$ is an integration constant. The dependence of
the canonical differential equations on $x$ and $y$ turns out to be also in
the most suitable form to reduce the repeated integrals of Eq. (56) to well-
known transcendental functions. In fact, the total differential of
$\mathbb{A}_{j}(x,y)$ is a $d$log-form of the type

$d\mathbb{A}=\sum_{i=1}^{12}\mathbb{M}_{i}\,d\log(\eta_{i})~{},$ (57)

with

$\begin{split}&\eta_{1}=x~{},\qquad\qquad\;\,\eta_{2}=1+x~{},\qquad\;\;\;\,\eta_{3}=1-x~{},\\\
&\eta_{4}=1+x^{2}~{},\qquad\;\eta_{5}=y~{},\qquad\qquad\;\;\,\eta_{6}=1+y~{},\\\
&\eta_{7}=1-y~{},\qquad\;\;\;\eta_{8}=1+y^{2}~{},\qquad\;\;\eta_{9}=x+y~{},\\\
&\eta_{10}=x-y~{},\qquad\eta_{11}=1+xy~{},\qquad\eta_{12}=1-xy~{}.\end{split}$
(58)

Upon factorization over the complex numbers of $\eta_{i}$, it is then possible
to express $\mathbf{I}^{(k)}(x,y)$ in terms of GPLs,

$G\left(\vec{a}_{n};t\right)=G\left(a_{1},\vec{a}_{n{-}1};t\right)=\int_{0}^{t}\frac{dt^{\prime}}{t^{\prime}-a_{1}}\,G\left(\vec{a}_{n{-}1};t^{\prime}\right),\qquad
G\left(\vec{0}_{n};t\right)=\frac{1}{n!}\,\log^{n}\left(t\right)\,.$ (59)

The iterative integration of Eq. (51) by means of Eq. (56) leaves, at every
order in the $\epsilon$ expansion, an undetermined boundary constant for each
MI. In principle, these constants can be determined by exploiting the
knowledge of the analytic expression of the MIs in specific kinematic limits.
Unfortunately, for the most complicated integrals, this information is
generally not easy to obtain through direct integration techniques.
Nevertheless, the differential equations are shown to also encode enough
information to determine most of the integration constants and, therefore,
minimize the needed amount of external independent information. In most cases,
the differential equations exhibit singularities, not only at the physical
thresholds of the MIs, but also in kinematics points which are regular points
of Feynman integrals (and, for this reason, often referred to as pseudo-
thresholds). The requirement of the solution of the differential equations
being regular at such nonphysical pseudo-thresholds can be translated into
linear equations for the integration constants. In particular, the boundary
constants have been determined in the following way:

  * •

The tadpoles $\text{I}_{1,2}$ and the factorized integral $\text{I}_{6}$ are
obtained by direct integration and provided as an external input to the system
of differential equations:

$\begin{split}&\text{I}_{1}=1\,,\\\
&\text{I}_{2}(x,y,\epsilon)=\left(\frac{(1-x^{2})^{2}\,y^{2}}{(1-y^{2})^{2}\,x^{2}}\right)^{-\epsilon}\,,\\\
&\text{I}_{6}(x,\epsilon)=\left(\frac{(1-x^{2})^{2}}{x^{2}}\right)^{-\epsilon}\,\left(1-\frac{\pi^{2}}{6}\,\epsilon^{2}-2\,\zeta_{3}\,\epsilon^{3}-\frac{\pi^{4}}{40}\,\epsilon^{4}+\mathcal{O}(\epsilon)^{5}\right)~{}.\end{split}$
(60)

  * •

The boundary constants of the integrals
$\text{I}_{3,4,5,9,10,11,12,14,16,17,18,19}$ are determined by imposing their
regularity at the pseudo-threshold $t\to 0$. In particular, due to the
prefactors that appear in the definition of the canonical MIs of Eq. (52),
$\text{I}_{3,4,5,9,10,14,16,17}$ vanish in this limit.999The fact that
$\text{I}_{14}$ vanishes in the $t\to 0$ limit can be inferred from the
results of Ref. [126]. The same conclusion is inferred for
$\text{I}_{11,12,18,19}$, analyzing the differential equation in this limit.
$\text{I}_{10}$ results to be vanishing in the $t\to 0$ limit, due to the
faster convergence of the factorized (canonical) massive bubble, with respect
to the massless one.

  * •

The boundary constants of $\text{I}_{7,8}$ are determined thanks to the
regularity at $m^{2}_{N}\to 0$; specifically, due to the prefactors in Eq.
(52), they vanish in this limit.

  * •

The boundary constants of $\text{I}_{13,15,20}$ are determined by demanding
regularity in the $t\to 4m^{2}_{N}$ limit. Thanks to the prefactors in Eq.
(52), they vanish in this limit.

The solution of Eq. (51) is originally derived in the region $0<x<1$ $\wedge$
$0<y<x$. In such region, all the polynomials $\eta_{i}$ are real and positive
and yield real-valued expression of the MIs. The definition of the kinematic
variables given in Eq. (15) can be inverted by choosing the solutions

$x=\frac{1}{2}\left(\sqrt{4-\sigma_{\ell}}-\sqrt{-\sigma_{\ell}}\right)~{},\qquad
y=\frac{1}{2}\left(\sqrt{4-\sigma_{N}}-\sqrt{-\sigma_{N}}\right)~{},\qquad\sigma_{i}=\frac{t}{m^{2}_{i}}~{},\qquad
i=N,\ell~{},$ (61)

so that the region $0<x<1$ $\wedge$ $0<y<x$ corresponds to the kinematic
configuration $t<0$ $\wedge$ $0<m^{2}_{N}<m^{2}_{\ell}$. At negative $t$, the
same definition of $x$ and $y$ is also valid in the reversed mass hierarchy
$0<m^{2}_{\ell}<m^{2}_{N}$ and therefore, can be used over the whole kinematic
region of interest for our the study of the form factors. However, once the
analytical expressions of the MIs are known in one region, analytic
continuation can be used to obtain the values of the MIs in any other
kinematic configuration of physical interest. In practice, this can be
achieved by propagating the Feynman prescription for positive values of the
square momentum transfer, $t\to t+i0^{+}$ to the $x$ and $y$ variables, so
that they acquire the proper imaginary part. All the GPLs expressions in this
work are manipulated via PolyLogTools [106] and an _in-house_ code. They are
numerically evaluated using GiNaC [105]. In some cases, high precision
numerical values for the boundary constants are fitted to analytic expressions
thanks to the PSLQ algorithm [127].

###  A.2 Master integrals for the equal-mass case

In this case ($m_{N}=m_{\ell}\equiv m$), all the scalar integrals are reduced
via IBPs to linear combinations of $15$ MIs. They obey a system of
differential equations which reads

$d\,\mathbf{I}(\epsilon,w)=\epsilon\,d\mathbb{A}(w)\,\mathbf{I}(\epsilon,w)\qquad\text{with}\qquad
df=dw\,\frac{\partial f}{\partial w}~{}.$ (62)

where the variable $w$ is related to the only scales present in the problem,
namely $t$ and $m^{2}$, via the relation

$-\frac{t}{m^{2}}=\frac{(1-w)^{2}}{w}~{}.$ (63)

A particular choice of MIs, for which the factorization in Eq. (62) holds, is
given by

$\displaystyle\text{I}_{1}$
$\displaystyle=\epsilon^{2}\,\mathcal{T}_{1}~{},\qquad$
$\displaystyle\text{I}_{2}=\epsilon^{2}\,\lambda_{m}\,\mathcal{T}_{2}~{},$
$\displaystyle\text{I}_{3}$
$\displaystyle=-\epsilon^{2}\,t\,\mathcal{T}_{3}~{},\qquad$
$\displaystyle\text{I}_{4}=\epsilon^{2}\,\left(\frac{1}{2}\,\mathcal{T}_{3}\,\left(\lambda_{m}+t\right)+\mathcal{T}_{4}\,\lambda_{m}\right)~{},$
$\displaystyle\text{I}_{5}$
$\displaystyle=-\epsilon^{2}\,t\,\mathcal{T}_{5}~{},\qquad$
$\displaystyle\text{I}_{6}=\epsilon^{2}\,m^{2}\,\mathcal{T}_{6}~{},$
$\displaystyle\text{I}_{7}$
$\displaystyle=-\epsilon^{2}\,t\,\lambda_{m}\,\mathcal{T}_{7}~{},\qquad$
$\displaystyle\text{I}_{8}=\epsilon^{3}\,\lambda_{m}\,\mathcal{T}_{8}~{},$
$\displaystyle\text{I}_{9}$
$\displaystyle=\epsilon^{2}\,\left(\frac{1}{4}\,\left(4\,m^{2}-\lambda_{m}-t\right)\,\left(\mathcal{T}_{3}+2\,\mathcal{T}_{4}\right)+m^{2}\,\left(4\,m^{2}-t\right)\,\mathcal{T}_{9}\right)~{},\qquad$
$\displaystyle\text{I}_{10}=\epsilon^{3}\,\lambda_{m}\,\mathcal{T}_{10}~{},$
$\displaystyle\text{I}_{11}$
$\displaystyle=\epsilon^{3}\,(1-2\,\epsilon)\,t\,\mathcal{T}_{11}~{},\qquad$
$\displaystyle\text{I}_{12}=\epsilon^{3}\,t\,\left(t-4\,m^{2}\right)\,\mathcal{T}_{12}~{},$
$\displaystyle\text{I}_{13}$
$\displaystyle=\epsilon^{3}\,\lambda_{m}\,\mathcal{T}_{13}~{},$
$\displaystyle\text{I}_{14}$
$\displaystyle=\epsilon^{2}\,\Big{(}\left(4\,m^{2}-\lambda_{m}-t\right)\,\left(\mathcal{T}_{2}-\epsilon\,\mathcal{T}_{13}\right)+(2\,\epsilon-1)\,\left(4\,m^{2}-t\right)\,\mathcal{T}_{14}\Big{)}~{},$
$\displaystyle\text{I}_{15}$
$\displaystyle=-\epsilon^{4}\,\lambda_{m}\,t\,\mathcal{T}_{15}~{},$ (64)

where we have introduced

$\lambda_{m}=\sqrt{-t}\,\sqrt{4\,m^{2}-t}~{}.$ (65)

(1)  $\mathcal{T}_{1}$ (2)  $\mathcal{T}_{2}$ (3)  $\mathcal{T}_{3}$ (4)
$\mathcal{T}_{4}$

(5)  $\mathcal{T}_{5}$ (6)  $\mathcal{T}_{6}$

(7)  $\mathcal{T}_{7}$ (8)  $\mathcal{T}_{8}$ (9)  $\mathcal{T}_{9}$ (10)
$\mathcal{T}_{10}$

(11)  $\mathcal{T}_{11}$ (12)  $\mathcal{T}_{12}$

(13)  $\mathcal{T}_{13}$ (14)  $\mathcal{T}_{14}$ (15)  $\mathcal{T}_{15}$

Figure 6: Master integrals required for the evaluation of the scalar and
pseudo-scalar form factors, corresponding to the two-loop diagrams in Fig. 1,
in the equal-mass case. In each diagram, we indicate the off-shell external
mediator (dashed lines), massive fermions with mass $m$ (blue lines), massless
photon (black lines) and squared denominators (dots).

The integral measure is defined as in Eq. (55), with the replacement
$m^{2}_{\ell}\to m^{2}$. The Taylor coefficients of the MIs,
$\mathbf{I}^{(k)}(w)$, can be obtained as described in Eqs. (54) and (56). In
this case, the total differential $d\mathbb{A}(w)$ reads

$d\,\mathbb{A}=\sum_{i=1}^{3}\mathbb{M}_{i}\,d\,\log\left(\eta_{i}\right)~{},$
(66)

where, once again, $\mathbb{M}_{i}$ are constant matrices and $n_{i}$ reduce
to

$\eta_{1}=w~{},\qquad\eta_{2}=1+w~{},\qquad\eta_{3}=1-w~{}.$ (67)

Thus, the iterated integrals introduced in Eq. (59) are written in terms of
$a_{i}=\pm 1,0$, yielding a more restricted class of functions, the _Harmonic
PolyLogarithms_ (HPLs) [102]. As in the general case, the fixing of the
boundary terms in Eq. (56) deserves a few comments.

  * •

The integrals $\text{I}_{1,5,6}$ are provided as an external input. In
particular, $\text{I}_{1,5}$ are obtained by direct integration

$\begin{split}&\text{I}_{1}=1~{},\\\
&\text{I}_{5}(w,\epsilon)=\left(\frac{(1-w)^{2}}{w}\right)^{-\epsilon}\,\left(1-\frac{\pi^{2}}{6}\,\epsilon^{2}-2\,\zeta_{3}\,\epsilon^{3}-\frac{\pi^{4}}{40}\,\epsilon^{4}+\mathcal{O}(\epsilon)^{5}\right)~{},\end{split}$
(68)

while $\text{I}_{6}$ is obtained from the general case as [128]

$\text{I}_{6}(\epsilon)=-\frac{\pi^{2}}{12}\,\epsilon^{2}+\frac{1}{4}\,\left(2\,\pi^{2}\,\log(2)-7\,\zeta_{3}\right)\,\epsilon^{3}+\frac{1}{360}\,\left(31\,\pi^{4}-180\,\log^{4}(2)-360\,\pi^{2}\,\log^{2}(2)-4320\,\text{Li}_{4}\left(\frac{1}{2}\right)\right)\,\epsilon^{4}+\mathcal{O}(\epsilon^{5})~{}.$
(69)

  * •

The integral $\text{I}_{12}$ is obtained as the product of $\text{I}_{2}$ and
$\text{I}_{10}$.

  * •

The boundary constants for $\text{I}_{2,3,4,8,9,11,13,14}$ are determined by
the regularity at the pseudo-threshold $t\to 0$. In particular
$\text{I}_{2,3,4,8,11,13}$ vanish in this limit due to the prefactors in Eq.
(64). Analyzing the differential equation in the $t\to 0$ limit, we obtain
relations among MIs in this limit, namely $\text{I}_{9}\left.\right|_{t\to
0}=3\,\text{I}_{6}$ and $\text{I}_{14}\left.\right|_{t\to 0}=6\,\text{I}_{6}$;
the latter are sufficient to determine the boundary constants for
$\text{I}_{9,14}$. $\text{I}_{7}$ vanishes due to the factorized canonical
bubble.

  * •

The boundary constants for $\text{I}_{10}$ are determined from the regularity
at the pseudo-threshold $t\to 4\,m^{2}$. In particular due to the prefactor in
Eq. (64), $\text{I}_{10}$ vanishes in this limit.

  * •

The boundary constants for $\text{I}_{10}$ are obtained by comparing our
results with those in Ref. [46].

Finally, we stress that Eq. (62) is solved in the region $0<w<1$, which
implies $t<0\,\wedge\,m^{2}>0$, where $\eta_{i}$’s are positive and the MIs
are real-valued. In such a regime, Eq. (63) can be inverted, yielding

$w=\frac{\sqrt{4\,m^{2}-t}-\sqrt{-t}}{\sqrt{4\,m^{2}-t}+\sqrt{-t}}~{}.$ (70)

##  Appendix B Limits of the scalar form factor,
$\bm{\mathcal{F}_{\mathcal{S}}^{\rm 1b}(t;m_{N},m_{\ell})}$

The expressions of the form factors for arbitrary kinematics are too long to
be shown here. Nevertheless, there are two limits which render manageable
analytical expressions. We hereby present the results for these two limiting
behaviors: the soft (i.e., $t\to 0$) and the equal-mass ($m_{N}\to m_{\ell}$)
limits.

###  B.1 Soft limit: $\bm{\mathcal{F}_{\mathcal{S}}^{\rm
1b}(t=0;m_{N},m_{\ell})}$

Figure 7: Diagrammatic representation of the $\bm{q^{\mu}\to 0}$ limit for the
one-body scalar form factor, $\bm{{\cal F}_{S}^{\rm 1b}}$. At zero momentum
transfer, no momentum flows into the diagram through the dashed line, which
just corresponds to the insertion of $\Lambda_{S}=\mathbb{I}$ in the fermionic
trace.

Setting $q^{\mu}=0$ before the integration of the loop momenta has the
advantage of drastically simplifying the whole IBPs reduction to MIs, as it
removes the dependence on $t$ from very the beginning. In order to understand
how the limit at the integrand level simplifies the reduction procedure, we
observe that the limit $q^{\mu}\to 0$ of Eq. (7) leaves only four independent
inverse propagators $D_{i}$ out of the six of the two-loop topology, since
$D_{2}\to D_{1}$ and $D_{5}\to D_{4}$. This implies that, at $q^{\mu}=0$, the
decomposition of the form factor into Feynman integrals given in Eq. (14)
reduces to

${\cal F}_{S}^{\rm
1b}\left(t=0;m_{N},m_{\ell}\right)=\sum_{\vec{n}}c_{\vec{n}\,j}(m_{i}^{2},\epsilon)\int\frac{{\differential^{d}k_{1}}}{(2\pi)^{d}}\frac{{\differential^{d}k_{2}}}{(2\pi)^{d}}\,\frac{D_{7}^{n_{7}}}{D_{1}^{n_{1}}D_{3}^{n_{3}}D_{4}^{n_{4}}D_{6}^{n_{6}}}\,,\quad
n_{i}\in\mathbb{N}~{},$ (71)

with

$\displaystyle D_{1}=k_{1}^{2}-m_{\ell}^{2}~{},\quad
D_{3}=k_{2}^{2}-m_{N}^{2}~{},\quad D_{4}=(k_{2}+p)^{2}~{},\quad
D_{6}=(k_{1}+k_{2}+p)^{2}-m_{\ell}^{2}~{},$ (72)

and where, as for arbitrary $t$, the auxiliary denominator
$D_{7}=(k_{1}-p)^{2}$ has been introduced so that all the five possible scalar
products between the loop momenta $k_{1,2}$ and the external momentum $p$ can
be expressed as linear combinations of $D_{i}$. Therefore, the $q^{\mu}\to 0$
limit can be understood, at a diagrammatic level, as the reduction of the two-
loop vertex function to the self-energy-like diagram depicted in Fig. 7. For
this case we identify $4$ MIs, which can be chosen as

$\displaystyle\text{I}_{1}=$
$\displaystyle\epsilon^{2}\,\mathcal{T}_{1}~{},\hskip 56.9055pt$
$\displaystyle\text{I}_{2}=\epsilon^{2}\,\mathcal{T}_{2}~{},$
$\displaystyle\text{I}_{3}=$
$\displaystyle\epsilon^{2}\,m_{\ell}\,m_{N}\left(\mathcal{T}_{3}+2\,\mathcal{T}_{4}\right)~{},\hskip
56.9055pt$
$\displaystyle\text{I}_{4}=\epsilon^{2}\,m^{2}_{N}\mathcal{T}_{4}~{},$ (73)

where $\mathcal{T}_{i}$’s are depicted in Fig. 8. These MIs obey a canonical
system of differential equations with respect to the variable $z_{N\ell}\equiv
m_{N}/m_{\ell}$, which can be solved in terms of HPLs. Normalizing the MIs as
in Eq. (55), the boundary conditions are obtained in the following way:

  * •

$\text{I}_{1}$ and $\text{I}_{2}$ are obtained by direct integration

$\displaystyle\text{I}_{1}$ $\displaystyle=$ $\displaystyle 1~{},$
$\displaystyle\text{I}_{2}(z,\epsilon)$ $\displaystyle=$ $\displaystyle
z^{-2\,\epsilon}~{}.$

  * •

The boundary conditions for $\text{I}_{3}$ and $\text{I}_{4}$ are determined
from to their regularity at $m_{N}\to 0$ [128]. In particular, thanks to the
prefactors in Eq. (74), they vanish in this limit.

The form factor ${\cal F}_{S}^{\rm 1b}(t=0;m_{N},m_{\ell})$ is then reduced
via IBPs to the following linear combination of MIs,

$\displaystyle{\cal F}_{S}^{\rm
1b}\left(t=0;m_{N},m_{\ell}\right)=-\frac{128\,\pi^{2}\,\alpha_{\text{em}}^{2}\,m_{\ell}}{\epsilon^{2}\,\left(\epsilon+1\right)\,\left(1-2\,\epsilon\right)\,\left(1-3\,\epsilon\right)\,m_{N}}\,\Big{[}$
$\displaystyle\epsilon\,\left(\epsilon\,\left(2\,\epsilon-1\right)+1\right)\,\text{I}_{1}+\epsilon\,\left(\epsilon\,\left(2\,\epsilon-3\right)-1\right)\,\text{I}_{2}$
$\displaystyle+\left(\frac{m_{N}}{m_{\ell}}\,\left(2\,\epsilon+1\right)\,\left(3\,\epsilon-1\right)-\frac{m_{\ell}}{m_{N}}\,\left(4\,\epsilon^{3}-4\,\epsilon^{2}-\epsilon+3\right)\right)\,\text{I}_{3}$
$\displaystyle+\left(\frac{m_{\ell}^{2}}{m_{N}^{2}}\,4\,\left(\epsilon-1\right)\,\left(\epsilon+1\right)\,\left(2\,\epsilon-1\right)-2\,\epsilon\,\left(\epsilon\,\left(2\,\epsilon+3\right)-3\right)\right)\,\text{I}_{4}\Big{]}~{}.$
(74)

(1)  $\mathcal{T}_{1}$ (2)  $\mathcal{T}_{2}$ (3)  $\mathcal{T}_{3}$ (4)
$\mathcal{T}_{4}$

Figure 8:  Master integrals required for the evaluation of the scalar form
factors in the soft limit, $\bm{t\to 0}$. The diagrams contain fermions with
mass $m_{\ell}$ (blue lines) and $m_{N}$ (red lines). External fermions are
on-shell.

After the insertion of the analytic expression of the four MIs, rewritten as
functions of the mass ratio $z_{N\ell}\equiv m_{N}/m_{\ell}$ only, restoring
the original integral measure and expanding Eq. (74) around four dimensions,
we obtain

$\displaystyle{\cal F}_{S}^{\rm
1b}\left(t=0;z_{N\ell}\equiv\frac{m_{N}}{m_{\ell}}\right)$ $\displaystyle=$
$\displaystyle-\frac{2\,\alpha_{\rm
em}^{2}}{\pi^{2}\,z_{N\ell}}\Big{\\{}1-\frac{\ln(z_{N\ell})}{2}+f_{S}(z_{N\ell})+f_{S}(-z_{N\ell})\Big{\\}}~{},$
(75) $\displaystyle\textrm{with}\hskip 19.91692ptf_{S}(\tau)$
$\displaystyle\equiv$
$\displaystyle\frac{1}{4\,\tau^{2}}\left(4+3\,\tau+\tau^{3}\right)\Big{[}\ln|\tau|\ln(1+\tau)+{\rm
Li}_{2}(-\tau)\Big{]}~{},$ (76)

as reported in Eq. (17).

###  B.2 Equal-mass limit: $\bm{\mathcal{F}_{\mathcal{S}}^{\rm 1b}(t;m,m)}$

In this limit, the results are expressed in terms of GPLs, $G(\vec{a},w)\equiv
G_{\vec{a}}(w)$, where $w$ was introduced in Eq. (63),

$\displaystyle\frac{\pi^{2}}{\alpha_{\rm em}^{2}}\,\mathcal{F}^{\rm
1b}_{\mathcal{S}}(w)=$
$\displaystyle-\frac{(20\,w^{2}+4\,w)}{\,(w-1)\,(w+1)^{2}}\,G_{0}(w)\;+\;\frac{w}{\,(w-1)\,(w+1)}\,\big{[}24\,G_{-1,0}(w)+4\,G_{1,0}(w)\big{]}$
$\displaystyle+\frac{w}{\,(w-1)^{2}\,(w+1)^{3}}[-12\,w^{3}-38\,w^{2}+8\,w+2]\,G_{0,0}(w)\;-\;\frac{w}{\pi^{2}\,(w-1)^{2}}\,\big{[}2\,G_{0,0,0}(w)+4\,\zeta_{3}\big{]}$
$\displaystyle+\frac{w\,\zeta_{3}}{\,(w-1)^{3}\,(w+1)}\,[5\,w^{2}-6\,w+5]\,G_{0}(w)\;+\;\frac{w}{\pi^{2}\,(w+1)^{2}}\,\left[12\,G_{0,0,1}(w)-\frac{4\,\pi^{2}}{3}\,G_{0}(w)+16\,G_{1}(w)\right]$
$\displaystyle+\frac{w}{\,(w-1)^{2}\,(w+1)^{2}}\,\Big{[}(-16\,w^{2}+64\,w-16)\,G_{0,-1,0}(w)+(8\,w^{2}-48\,w+8)\,G_{1,0,0}(w)+16\,w\,G_{0,1,0}(w)\Big{]}$
$\displaystyle+\frac{w}{\,(w-1)\,(w+1)^{3}}\Bigg{[}(-14\,w^{2}+4\,w-14)\,G_{0,1}(w)+\left(\frac{w^{2}}{3}+2\,w+\frac{1}{3}\right)\,G_{1,0}(w)\,\pi^{2}$
$\displaystyle+(2\,w^{2}+12\,w+2)\,\big{(}G_{1,0,0,0}(w)\,-\,G_{0,0,1,0}(w)\big{)}+\left(2\,w^{2}-\frac{20\,w}{3}+2\right)\pi^{2}\Bigg{]}$
$\displaystyle+\frac{w}{\,(w-1)^{3}\,(w+1)^{3}}\Bigg{[}\left(\frac{5\,w^{4}}{6}-2\,w^{3}+\frac{23\,w^{2}}{3}-2\,w+\frac{5}{6}\right)\,G_{0,0}(w)\,\pi^{2}$
$\displaystyle+\left(8\,w^{4}-32\,w^{3}+112\,w^{2}-32\,w+8\right)\,G_{0,0,-1,0}(w)+\left(\frac{3\,w^{4}}{2}-2\,w^{3}+9\,w^{2}-2\,w+\frac{3}{2}\right)\,G_{0,0,0,0}(w)$
$\displaystyle+\left(-9\,w^{4}+12\,w^{3}-54\,w^{2}+12\,w-9\right)\,G_{0,0,0,1}(w)+\left(-2\,w^{4}+24\,w^{3}-76\,w^{2}+24\,w-2\right)\,G_{0,1,0,0}(w)$
$\displaystyle+\left(\frac{w^{4}}{18}-\frac{2\,w^{3}}{45}+\frac{11\,w^{2}}{45}-\frac{2\,w}{45}+\frac{1}{18}\right)\,\pi^{4}\Bigg{]}~{}.$
(77)

This form factor is related, by crossing symmetry, to the one appearing in the
amplitude for a (scalar) Higgs boson decaying into a heavy-quark pair [14],
which contains the expression of the total form factor, due to the
contribution of several vertex diagrams, in $s$-channel. The expression given
above is the contribution of a single diagram, namely the vertex diagram with
a closed triangle-loop of a massive fermion, in $t$-channel, and it is shown
here, for the first time.

##  Appendix C Limits of the pseudo-scalar form factor,
$\bm{\mathcal{F}_{\mathcal{P}}^{\rm 1b}(t;m_{N},m_{\ell})}$

As for the scalar case, we provide here the expressions for the soft (i.e.,
$t\sim 0$) and the equal-mass (i.e., $m_{N}\to m_{\ell}$) limits for the
pseudo-scalar form factor.

###  C.1 Soft asymptotic limit: $\bm{\mathcal{F}_{\mathcal{P}}^{\rm 1b}(t\sim
0;m_{N},m_{\ell})}$

We consider an asymptotic expansion around $t\sim 0$ for the full integrated
result, by exploiting the algebraic properties and functional relations among
GPLs. The results are expressed in terms of the dimensionless variables
$(x,y)$ introduced in Eq. (15). Notice that $t\sim 0$ and fixed $m_{N,\ell}$
corresponds to $(x,y)\sim(1,1)$. The ordering $x<y$ or $y<x$ is kept fixed in
this limit, and depends on the fermion mass hierarchy.

For the mass hierarchy $m_{l}^{2}<m_{N}^{2}$ (i.e., $0<x<y$):

$\displaystyle\frac{\pi^{2}}{\alpha_{\rm
em}^{2}}\,\mathcal{F}_{\mathcal{P}}^{\rm 1b}(t\sim 0;m_{N}>m_{\ell})=$
$\displaystyle-\frac{\pi^{2}}{3}+\frac{(1-x)}{\,(1-y)}\,\big{[}-2\,\log(1-x)-2\,\log(2)+3\big{]}+\frac{(1-x)^{2}}{\,(1-y)}\,\left[-\log(1-x)-\log(2)+\frac{1}{2}\right]$
$\displaystyle+(1-x)^{2}\,\left[\frac{\log(1-x)}{2}+\frac{\log(2)}{2}+\frac{37\,\pi^{2}}{360}-\frac{1}{4}\right]+\frac{(1-x)}{\pi^{2}}\,\left[\log(1-x)+\log(2)+\frac{\pi^{2}}{4}-\frac{3}{2}\right]$
$\displaystyle+\mathcal{O}\big{(}(1-y),\;(1-x)^{3}\big{)}~{}.$ (78)

Here, $\mathcal{O}\big{(}(1-y)^{n},\;(1-x)^{m}\big{)}$ indicates that the
expression was first expanded in $y$ to order $n-1$, and then expanded in $x$
to order $m-1$.

For the mass hierarchy $m_{N}^{2}<m_{l}^{2}$ (i.e., $0<y<x$):

$\displaystyle\frac{\pi^{2}}{\alpha_{\rm
em}^{2}}\,\mathcal{F}_{\mathcal{P}}^{\rm 1b}(t\sim 0;m_{N}<m_{\ell})=$
$\displaystyle\frac{(1-x)}{\,(1-y)}\,\left[-\frac{3\log(1-x)}{2}-\frac{\log(1-y)}{2}-2\,\log(2)+\frac{17}{8}\right]$
$\displaystyle+(1-x)\,\left[\frac{3\,\log(1-x)}{4}+\frac{\log(1-y)}{4}+\log(2)+\frac{\pi^{2}}{4}-\frac{21}{16}\right]+\mathcal{O}\big{(}(1-x)^{2},\;(1-y)\big{)}~{}.$
(79)

Analogously, $\mathcal{O}\big{(}(1-x)^{n},\;(1-y)^{m}\big{)}$ indicates that
the expression was first expanded in $x$ to order $n-1$, and then expanded in
$y$ to order $m-1$.

###  C.2 Equal-mass limit: $\bm{\mathcal{F}_{\mathcal{P}}^{\rm 1b}(t;m,m)}$

The results for the equal-mass limit are expressed in terms GPLs:
$G(\vec{a},w)\equiv G_{\vec{a}}(w)$, where $w$ was introduced in Eq. (63),

$\displaystyle\frac{\pi^{2}}{\alpha_{\rm em}^{2}}\,\mathcal{F}^{\rm
1b}_{\mathcal{P}}(w)=$
$\displaystyle-\frac{w^{2}}{\,(w-1)^{3}\,(w+1)}\,\left[\frac{2\pi^{2}}{3}\,G_{0}(w)+4\,G_{0,0,0}(w)\right]$
$\displaystyle+\frac{w}{\,(w-1)^{2}}\,\left[G_{0,0}(w)-4\,G_{0,0,1}(w)-4\,G_{0,1,0}(w)+8\,G_{1,0,0}(w)+12\,\zeta_{3}+\frac{\pi^{2}}{3}\right]$
$\displaystyle+\frac{w}{\,(w-1)\,(w+1)}\,\left[-3\,\zeta_{3}\,G_{0}(w)-\frac{\pi^{2}}{6}\,G_{0,0}(w)-\frac{\pi^{2}}{3}\,G_{1,0}(w)-\frac{1}{2}\,G_{0,0,0,0}(w)\right.$
$\displaystyle\hskip
89.626pt\left.+3\,G_{0,0,0,1}(w)+2\,G_{0,0,1,0}(w)-2\,G_{0,1,0,0}(w)-2\,G_{1,0,0,0}(w)-\frac{\pi^{4}}{45}\right]~{}.$
(80)

The calculation of the two-loop pseudo-scalar form factor in the $s$-channel
was considered previously, for the two-loop QCD corrections to the heavy quark
form factors (see the expression of $A_{R}$ in Ref. [13]). We have reproduced
it using our MIs, which are, instead, defined for a $t$-channel process. The
result shown above refers to a $t$-channel process, and it is presented here
for the first time.

## References

  * Argeri and Mastrolia [2007] M. Argeri and P. Mastrolia, Feynman diagrams and differential equations, Int. J. Mod. Phys. A22, 4375 (2007), arXiv:0707.4037 [hep-ph] .
  * Grozin [2011] A. G. Grozin, Integration by parts: an introduction, Int. J. Mod. Phys. A26, 2807 (2011), arXiv:1104.3993 [hep-ph] .
  * Smirnov [2012] V. A. Smirnov, Analytic tools for Feynman integrals, Springer Tracts Mod. Phys. 250, 1 (2012).
  * Henn [2015] J. M. Henn, Lectures on differential equations for Feynman integrals, J. Phys. A48, 153001 (2015), arXiv:1412.2296 [hep-ph] .
  * Duhr [2014] C. Duhr, Mathematical aspects of scattering amplitudes, in  _Theoretical Advanced Study Institute in Elementary Particle Physics: Journeys Through the Precision Frontier: Amplitudes for Colliders (TASI 2014) Boulder, Colorado, June 2-27, 2014_ (2014) arXiv:1411.7538 [hep-ph] .
  * Altarelli and Lampe [1993] G. Altarelli and B. Lampe, Second order QCD corrections to heavy quark forward-backward asymmetries, Nucl. Phys. B391, 3 (1993).
  * Ravindran and van Neerven [1998] V. Ravindran and W. L. van Neerven, Second order QCD corrections to the forward-backward asymmetry in e+ e- collisions, Phys. Lett. B445, 214 (1998), arXiv:hep-ph/9809411 [hep-ph] .
  * Catani and Seymour [1999] S. Catani and M. H. Seymour, Corrections of $\mathcal{O}(\alpha_{s}^{2})$ to the forward backward asymmetry, JHEP 07, 023, arXiv:hep-ph/9905424 [hep-ph] .
  * Bernreuther  _et al._ [2005a] W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, T. Leineweber, P. Mastrolia, and E. Remiddi, Two-loop QCD corrections to the heavy quark form-factors: the vector contributions, Nucl. Phys. B706, 245 (2005a), arXiv:hep-ph/0406046 [hep-ph] .
  * Bernreuther  _et al._ [2005b] W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, T. Leineweber, P. Mastrolia, and E. Remiddi, Two-loop QCD corrections to the heavy quark form-factors: axial vector contributions, Nucl. Phys. B712, 229 (2005b), arXiv:hep-ph/0412259 [hep-ph] .
  * Bonciani  _et al._ [2008] R. Bonciani, A. Ferroglia, T. Gehrmann, D. Maitre, and C. Studerus, Two-loop fermionic corrections to heavy-quark pair production: the quark-antiquark channel, JHEP 07, 129, arXiv:0806.2301 [hep-ph] .
  * Ablinger  _et al._ [2018] J. Ablinger, A. Behring, J. Blümlein, G. Falcioni, A. De Freitas, P. Marquard, N. Rana, and C. Schneider, Heavy quark form factors at two loops, Phys. Rev. D97, 094022 (2018), arXiv:1712.09889 [hep-ph] .
  * Bernreuther  _et al._ [2005c] W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, T. Leineweber, and E. Remiddi, Two-loop QCD corrections to the heavy quark form-factors: anomaly contributions, Nucl. Phys. B723, 91 (2005c), arXiv:hep-ph/0504190 [hep-ph] .
  * Bernreuther  _et al._ [2005d] W. Bernreuther, R. Bonciani, T. Gehrmann, R. Heinesch, P. Mastrolia, and E. Remiddi, Decays of scalar and pseudoscalar Higgs bosons into fermions: two-loop QCD corrections to the Higgs-quark-antiquark amplitude, Phys. Rev. D72, 096002 (2005d), arXiv:hep-ph/0508254 [hep-ph] .
  * Bernreuther  _et al._ [2018] W. Bernreuther, L. Chen, and Z.-G. Si, Differential decay rates of CP-even and CP-odd Higgs bosons to top and bottom quarks at NNLO QCD, JHEP 07, 159, arXiv:1805.06658 [hep-ph] .
  * Primo  _et al._ [2019] A. Primo, G. Sasso, G. Somogyi, and F. Tramontano, Exact top Yukawa corrections to Higgs boson decay into bottom quarks, Phys. Rev. D99, 054013 (2019), arXiv:1812.07811 [hep-ph] .
  * Behring and Bizoń [2020] A. Behring and W. Bizoń, Higgs decay into massive b-quarks at NNLO QCD in the nested soft-collinear subtraction scheme, JHEP 01, 189, arXiv:1911.11524 [hep-ph] .
  * Mondini  _et al._ [2020] R. Mondini, U. Schubert, and C. Williams, Top-induced contributions to $H\rightarrow b\bar{b}$ and $H\rightarrow c\bar{c}$ at $\mathcal{O}(\alpha_{s}^{3})$, JHEP 12, 058, arXiv:2006.03563 [hep-ph] .
  * Heinemeyer  _et al._ [2004] S. Heinemeyer, D. Stockinger, and G. Weiglein, Electroweak and supersymmetric two-loop corrections to (g-2)(mu), Nucl. Phys. B699, 103 (2004), arXiv:hep-ph/0405255 [hep-ph] .
  * Stockinger [2007] D. Stockinger, The muon magnetic moment and supersymmetry, J. Phys. G34, R45 (2007), arXiv:hep-ph/0609168 [hep-ph] .
  * Miller  _et al._ [2007] J. P. Miller, E. de Rafael, and B. L. Roberts, Muon (g-2): experiment and theory, Rept. Prog. Phys. 70, 795 (2007), arXiv:hep-ph/0703049 [hep-ph] .
  * Davidson and Grenier [2010] S. Davidson and G. J. Grenier, Lepton flavour violating Higgs and tau to mu gamma, Phys. Rev. D81, 095016 (2010), arXiv:1001.0434 [hep-ph] .
  * Ilisie [2015] V. Ilisie, New Barr-Zee contributions to $(g-2)_{\mu}$ in two-Higgs-doublet models, JHEP 04, 077, arXiv:1502.04199 [hep-ph] .
  * Hisano  _et al._ [2018] J. Hisano, R. Nagai, and N. Nagata, Singlet Dirac fermion dark matter with mediators at loop, JHEP 12, 059, arXiv:1808.06301 [hep-ph] .
  * Dong  _et al._ [2020] X.-X. Dong, S.-M. Zhao, H.-B. Zhang, and T.-F. Feng, The two-loop corrections to lepton MDMs and EDMs in the EBLMSSM, J. Phys. G 47, 045002 (2020), arXiv:1901.07701 [hep-ph] .
  * Barr and Zee [1990] S. M. Barr and A. Zee, Electric dipole moment of the electron and of the neutron, Phys. Rev. Lett. 65, 21 (1990), [Erratum: Phys. Rev. Lett. 65 (1990) 2920].
  * Caswell and Lepage [1986] W. E. Caswell and G. P. Lepage, Effective lagrangians for bound state problems in QED, QCD, and other field theories, Phys. Lett. 167B, 437 (1986).
  * Brambilla  _et al._ [2005] N. Brambilla, A. Pineda, J. Soto, and A. Vairo, Effective field theories for heavy quarkonium, Rev. Mod. Phys. 77, 1423 (2005), arXiv:hep-ph/0410047 [hep-ph] .
  * Pineda [2012] A. Pineda, Review of heavy quarkonium at weak coupling, Prog. Part. Nucl. Phys. 67, 735 (2012), arXiv:1111.0165 [hep-ph] .
  * Gerlach  _et al._ [2019] M. Gerlach, G. Mishima, and M. Steinhauser, Matching coefficients in nonrelativistic QCD to two-loop accuracy, Phys. Rev. D100, 054016 (2019), arXiv:1907.08227 [hep-ph] .
  * Pratap and Smith [1972] M. Pratap and J. Smith, Eta and pion decays into lepton pairs, Phys. Rev. D5, 2020 (1972).
  * Bergstrom  _et al._ [1983] L. Bergstrom, E. Masso, L. Ametller, and A. Bramon, $Q^{2}$ duality and rare pion decays, Phys. Lett. 126B, 117 (1983).
  * Ametller  _et al._ [1983] L. Ametller, L. Bergstrom, A. Bramon, and E. Masso, The quark triangle: application to pion and $\eta$ decays, Nucl. Phys. B228, 301 (1983).
  * Pich and Bernabeu [1984] A. Pich and J. Bernabeu, Rare decay modes of the neutral pion, Z. Phys. C22, 197 (1984).
  * Gorishnii  _et al._ [1990] S. G. Gorishnii, A. L. Kataev, S. A. Larin, and L. R. Surguladze, Corrected three Loop QCD correction to the correlator of the quark scalar currents and $\gamma$ (Tot) ($H^{0}\to$ hadrons), Mod. Phys. Lett. A5, 2703 (1990).
  * Surguladze [1994a] L. R. Surguladze, Quark mass effects in fermionic decays of the Higgs boson in ${\cal O}(\alpha_{s}^{2})$ perturbative QCD, Phys. Lett. B341, 60 (1994a), arXiv:hep-ph/9405325 [hep-ph] .
  * Surguladze [1994b] L. R. Surguladze, Minimal supersymmetric Higgs boson decay rate in ${\cal O}(\alpha_{s}^{2})$ perturbative QCD, Phys. Lett. B338, 229 (1994b), arXiv:hep-ph/9406294 [hep-ph] .
  * Chetyrkin and Kwiatkowski [1996] K. G. Chetyrkin and A. Kwiatkowski, Second order QCD corrections to scalar and pseudoscalar Higgs decays into massive bottom quarks, Nucl. Phys. B461, 3 (1996), arXiv:hep-ph/9505358 [hep-ph] .
  * Chetyrkin  _et al._ [1997] K. G. Chetyrkin, B. A. Kniehl, and M. Steinhauser, Three loop $\mathcal{O}(\alpha_{s}^{2}\,G_{F}\,M_{t}^{2})$ corrections to hadronic Higgs decays, Nucl. Phys. B490, 19 (1997), arXiv:hep-ph/9701277 [hep-ph] .
  * Larin  _et al._ [1995] S. A. Larin, T. van Ritbergen, and J. A. M. Vermaseren, The large top quark mass expansion for Higgs boson decays into bottom quarks and into gluons, Phys. Lett. B362, 134 (1995), arXiv:hep-ph/9506465 [hep-ph] .
  * Harlander and Steinhauser [1997] R. Harlander and M. Steinhauser, Higgs decay to top quarks at $\mathcal{O}(\alpha_{s}^{2})$, Phys. Rev. D56, 3980 (1997), arXiv:hep-ph/9704436 [hep-ph] .
  * Chetyrkin  _et al._ [1998] K. G. Chetyrkin, R. Harlander, and M. Steinhauser, Singlet polarization functions at $\mathcal{O}(\alpha_{s}^{2})$, Phys. Rev. D58, 014012 (1998), arXiv:hep-ph/9801432 [hep-ph] .
  * Kotikov [1991] A. V. Kotikov, Differential equations method: new technique for massive Feynman diagrams calculation, Phys. Lett. B254, 158 (1991).
  * Fleischer  _et al._ [1999] J. Fleischer, A. V. Kotikov, and O. L. Veretin, Analytic two loop results for selfenergy type and vertex type diagrams with one nonzero mass, Nucl. Phys. B547, 343 (1999), arXiv:hep-ph/9808242 [hep-ph] .
  * Gehrmann and Remiddi [2000] T. Gehrmann and E. Remiddi, Differential equations for two loop four point functions, Nucl. Phys. B580, 485 (2000), arXiv:hep-ph/9912329 [hep-ph] .
  * Bonciani  _et al._ [2004] R. Bonciani, P. Mastrolia, and E. Remiddi, Master integrals for the two loop QCD virtual corrections to the forward backward asymmetry, Nucl. Phys. B690, 138 (2004), arXiv:hep-ph/0311145 [hep-ph] .
  * Davydychev and Kalmykov [2004] A. I. Davydychev and M. Yu. Kalmykov, Massive Feynman diagrams and inverse binomial sums, Nucl. Phys. B699, 3 (2004), arXiv:hep-th/0303162 [hep-th] .
  * Bernabei  _et al._ [2008] R. Bernabei _et al._ , Investigating electron interacting dark matter, Phys. Rev. D77, 023506 (2008), arXiv:0712.0562 [astro-ph] .
  * Dedes  _et al._ [2010] A. Dedes, I. Giomataris, K. Suxho, and J. D. Vergados, Searching for secluded dark matter via direct detection of recoiling nuclei as well as low energy electrons, Nucl. Phys. B826, 148 (2010), arXiv:0907.0758 [hep-ph] .
  * Kopp  _et al._ [2009] J. Kopp, V. Niro, T. Schwetz, and J. Zupan, DAMA/LIBRA and leptonically interacting dark matter, Phys. Rev. D80, 083502 (2009), arXiv:0907.3159 [hep-ph] .
  * Feldstein  _et al._ [2010] B. Feldstein, P. W. Graham, and S. Rajendran, Luminous dark matter, Phys. Rev. D82, 075019 (2010), arXiv:1008.1988 [hep-ph] .
  * Chang  _et al._ [2014] S. Chang, R. Edezhath, J. Hutchinson, and M. Luty, Leptophilic effective WIMPs, Phys. Rev. D90, 015011 (2014), arXiv:1402.7358 [hep-ph] .
  * Bell  _et al._ [2014] N. F. Bell, Y. Cai, R. K. Leane, and A. D. Medina, Leptophilic dark matter with $Z^{\prime}$ interactions, Phys. Rev. D90, 035027 (2014), arXiv:1407.3001 [hep-ph] .
  * Foot [2014] R. Foot, Can dark matter-electron scattering explain the DAMA annual modulation signal?, Phys. Rev. D90, 121302 (2014), arXiv:1407.4213 [hep-ph] .
  * Roberts  _et al._ [2016] B. M. Roberts, V. A. Dzuba, V. V. Flambaum, M. Pospelov, and Y. V. Stadnik, Dark matter scattering on electrons: Accurate calculations of atomic excitations and implications for the DAMA signal, Phys. Rev. D93, 115037 (2016), arXiv:1604.04559 [hep-ph] .
  * Fox and Poppitz [2009] P. J. Fox and E. Poppitz, Leptophilic dark matter, Phys. Rev. D79, 083528 (2009), arXiv:0811.0399 [hep-ph] .
  * Cao  _et al._ [2009] Q.-H. Cao, E. Ma, and G. Shaughnessy, Dark matter: the leptonic connection, Phys. Lett. B673, 152 (2009), arXiv:0901.1334 [hep-ph] .
  * Bi  _et al._ [2009] X.-J. Bi, X.-G. He, and Q. Yuan, Parameters in a class of leptophilic models from PAMELA, ATIC and FERMI, Phys. Lett. B678, 168 (2009), arXiv:0903.0122 [hep-ph] .
  * Ibarra  _et al._ [2009] A. Ibarra, A. Ringwald, D. Tran, and C. Weniger, Cosmic rays from leptophilic dark matter decay via kinetic mixing, JCAP 0908, 017, arXiv:0903.3625 [hep-ph] .
  * Cohen and Zurek [2010] T. Cohen and K. M. Zurek, Leptophilic dark matter from the lepton asymmetry, Phys. Rev. Lett. 104, 101301 (2010), arXiv:0909.2035 [hep-ph] .
  * Cavasonza  _et al._ [2017] L. A. Cavasonza, H. Gast, M. Krämer, M. Pellen, and S. Schael, Constraints on leptophilic dark matter from the AMS-02 experiment, Astrophys. J. 839, 36 (2017), [Erratum: Astrophys. J. 869 (2018) 89], arXiv:1612.06634 [hep-ph] .
  * Duan  _et al._ [2018a] G. H. Duan, L. Feng, F. Wang, L. Wu, J. M. Yang, and R. Zheng, Simplified TeV leptophilic dark matter in light of DAMPE data, JHEP 02, 107, arXiv:1711.11012 [hep-ph] .
  * Athron  _et al._ [2018] P. Athron, C. Balazs, A. Fowlie, and Y. Zhang, Model-independent analysis of the DAMPE excess, JHEP 02, 121, arXiv:1711.11376 [hep-ph] .
  * Duan  _et al._ [2018b] G. H. Duan, X.-G. He, L. Wu, and J. M. Yang, Leptophilic dark matter in gauged $U(1)_{L{{}_{e}}-L_{\mu}}$ model in light of DAMPE cosmic ray ${e{{}^{+}}}+{e{{}^{-}}}$ excess, Eur. Phys. J. C78, 323 (2018b), arXiv:1711.11563 [hep-ph] .
  * Ghorbani and Ghorbani [2018] K. Ghorbani and P. H. Ghorbani, DAMPE electron-positron excess in leptophilic $Z^{\prime}$ model, JHEP 05, 125, arXiv:1712.01239 [hep-ph] .
  * Han  _et al._ [2018] Z.-L. Han, W. Wang, and R. Ding, Radiative seesaw model and DAMPE excess from leptophilic gauge symmetry, Eur. Phys. J. C78, 216 (2018), arXiv:1712.05722 [hep-ph] .
  * Yaser Ayazi and Mohamadnejad [2020] S. Yaser Ayazi and A. Mohamadnejad, DAMPE excess from leptophilic vector dark matter: a model-independent approach, J. Phys. G 47, 095003 (2020), arXiv:1909.10729 [hep-ph] .
  * Lu and Zong [2016] B.-Q. Lu and H.-S. Zong, Leptophilic dark matter in galactic center excess, Phys. Rev. D93, 083504 (2016), [Addendum: Phys. Rev. D93 (2016) 089910].
  * Agrawal  _et al._ [2014] P. Agrawal, Z. Chacko, and C. B. Verhaaren, Leptophilic dark matter and the anomalous magnetic moment of the muon, JHEP 08, 147, arXiv:1402.7369 [hep-ph] .
  * Bandyopadhyay  _et al._ [2018] P. Bandyopadhyay, E. J. Chun, and R. Mandal, Scalar dark matter in leptophilic two-Higgs-doublet model, Phys. Lett. B779, 201 (2018), arXiv:1709.08581 [hep-ph] .
  * Chen  _et al._ [2018] C.-Y. Chen, J. Kozaczuk, and Y.-M. Zhong, Exploring leptophilic dark matter with NA64-$\mu$, JHEP 10, 154, arXiv:1807.03790 [hep-ph] .
  * Buckley and Feld [2015] M. R. Buckley and D. Feld, Dark matter in leptophilic Higgs models after the LHC Run-I, Phys. Rev. D92, 075024 (2015), arXiv:1508.00908 [hep-ph] .
  * D’Eramo  _et al._ [2017] F. D’Eramo, B. J. Kavanagh, and P. Panci, Probing leptophilic dark sectors with hadronic processes, Phys. Lett. B771, 339 (2017), arXiv:1702.00016 [hep-ph] .
  * Dutta  _et al._ [2017] S. Dutta, D. Sachdeva, and B. Rawat, Signals of leptophilic dark matter at the ILC, Eur. Phys. J. C77, 639 (2017), arXiv:1704.03994 [hep-ph] .
  * Madge and Schwaller [2019] E. Madge and P. Schwaller, Leptophilic dark matter from gauged lepton number: Phenomenology and gravitational wave signatures, JHEP 02, 048, arXiv:1809.09110 [hep-ph] .
  * Marsicano  _et al._ [2018] L. Marsicano, M. Battaglieri, A. Celentano, R. De Vita, and Y.-M. Zhong, Probing leptophilic dark sectors at electron beam-dump facilities, Phys. Rev. D98, 115022 (2018), arXiv:1812.03829 [hep-ex] .
  * Junius  _et al._ [2019] S. Junius, L. Lopez-Honorez, and A. Mariotti, A feeble window on leptophilic dark matter, JHEP 07, 136, arXiv:1904.07513 [hep-ph] .
  * Garani and Palomares-Ruiz [2017] R. Garani and S. Palomares-Ruiz, Dark matter in the Sun: scattering off electrons vs nucleons, JCAP 1705 (05), 007, arXiv:1702.02768 [hep-ph] .
  * Liang  _et al._ [2018] Z.-L. Liang, Y.-L. Tang, and Z.-Q. Yang, The leptophilic dark matter in the Sun: the minimum testable mass, JCAP 1810 (10), 035, arXiv:1802.01005 [hep-ph] .
  * Bell  _et al._ [2019] N. F. Bell, G. Busoni, and S. Robles, Capture of leptophilic dark matter in neutron stars, JCAP 1906 (06), 054, arXiv:1904.09803 [hep-ph] .
  * Garani and Heeck [2019] R. Garani and J. Heeck, Dark matter interactions with muons in neutron stars, Phys. Rev. D100, 035039 (2019), arXiv:1906.10145 [hep-ph] .
  * Joglekar  _et al._ [2020a] A. Joglekar, N. Raj, P. Tanedo, and H.-B. Yu, Relativistic capture of dark matter by electrons in neutron stars, Phys. Lett. B809, 135767 (2020a), arXiv:1911.13293 [hep-ph] .
  * Joglekar  _et al._ [2020b] A. Joglekar, N. Raj, P. Tanedo, and H.-B. Yu, Dark kinetic heating of neutron stars from contact interactions with relativistic targets, Phys. Rev. D102, 123002 (2020b), arXiv:2004.09539 [hep-ph] .
  * Bell  _et al._ [2021] N. F. Bell, G. Busoni, S. Robles, and M. Virgato, Improved rreatment of dark matter capture in neutron stars II: leptonic targets, JCAP 03, 086, arXiv:2010.13257 [hep-ph] .
  * Essig  _et al._ [2012] R. Essig, J. Mardon, and T. Volansky, Direct detection of sub-GeV dark matter, Phys. Rev. D85, 076007 (2012), arXiv:1108.5383 [hep-ph] .
  * Chen  _et al._ [2015] J.-W. Chen, H.-C. Chi, C. P. Liu, C.-L. Wu, and C.-P. Wu, Electronic and nuclear contributions in sub-GeV dark matter scattering: A case study with hydrogen, Phys. Rev. D92, 096013 (2015), arXiv:1508.03508 [hep-ph] .
  * Lee  _et al._ [2015] S. K. Lee, M. Lisanti, S. Mishra-Sharma, and B. R. Safdi, Modulation effects in dark matter-electron scattering experiments, Phys. Rev. D92, 083517 (2015), arXiv:1508.07361 [hep-ph] .
  * Essig  _et al._ [2016] R. Essig _et al._ , Direct detection of sub-GeV dark matter with semiconductor targets, JHEP 05, 046, arXiv:1509.01598 [hep-ph] .
  * Cappiello  _et al._ [2019] C. V. Cappiello, K. C. Y. Ng, and J. F. Beacom, Reverse direct detection: cosmic ray scattering with light dark matter, Phys. Rev. D99, 063004 (2019), arXiv:1810.07705 [hep-ph] .
  * An  _et al._ [2018] H. An, M. Pospelov, J. Pradler, and A. Ritz, Directly detecting MeV-scale dark matter via solar reflection, Phys. Rev. Lett. 120, 141801 (2018), [Erratum: Phys. Rev. Lett. 121, no.25 (2018) 259903], arXiv:1708.03642 [hep-ph] .
  * Ema  _et al._ [2019] Y. Ema, F. Sala, and R. Sato, Light dark matter at neutrino experiments, Phys. Rev. Lett. 122, 181802 (2019), arXiv:1811.00520 [hep-ph] .
  * Cappiello and Beacom [2019] C. Cappiello and J. F. Beacom, Strong new limits on light dark matter from neutrino experiments, Phys. Rev. D100, 103011 (2019), arXiv:1906.11283 [hep-ph] .
  * Frandsen  _et al._ [2012] M. T. Frandsen, U. Haisch, F. Kahlhoefer, P. Mertsch, and K. Schmidt-Hoberg, Loop-induced dark matter direct detection signals from gamma-ray lines, JCAP 1210, 033, arXiv:1207.3971 [hep-ph] .
  * Weiner and Yavin [2012] N. Weiner and I. Yavin, How dark are Majorana WIMPs? Signals from MiDM and Rayleigh dark matter, Phys. Rev. D86, 075021 (2012), arXiv:1206.2910 [hep-ph] .
  * Weiner and Yavin [2013] N. Weiner and I. Yavin, UV completions of magnetic inelastic and Rayleigh dark matter for the Fermi line(s), Phys. Rev. D87, 023523 (2013), arXiv:1209.1093 [hep-ph] .
  * Kavanagh  _et al._ [2019] B. J. Kavanagh, P. Panci, and R. Ziegler, Faint light from dark matter: classifying and constraining dark matter-photon effective operators, JHEP 04, 089, arXiv:1810.00033 [hep-ph] .
  * Ovanesyan and Vecchi [2015] G. Ovanesyan and L. Vecchi, Direct detection of dark matter polarizability, JHEP 07, 128, arXiv:1410.0601 [hep-ph] .
  * Chetyrkin and Tkachov [1981] K. G. Chetyrkin and F. V. Tkachov, Integration by parts: the algorithm to calculate beta functions in 4 loops, Nucl. Phys. B192, 159 (1981).
  * Laporta [2000] S. Laporta, High precision calculation of multiloop Feynman integrals by difference equations, Int. J. Mod. Phys. A 15, 5087 (2000), arXiv:hep-ph/0102033 .
  * [100] A. von Manteuffel and C. Studerus, Reduze 2 - Distributed Feynman Integral Reduction, arXiv:1201.4330 [hep-ph] .
  * Goncharov [1995] A. Goncharov, Polylogarithms in arithmetic and geometry, Proceedings of the International Congree of Mathematicians 1,2, 374 (1995).
  * Remiddi and Vermaseren [2000] E. Remiddi and J. Vermaseren, Harmonic polylogarithms, Int. J. Mod. Phys. A15, 725 (2000), arXiv:hep-ph/9905237 .
  * Gehrmann and Remiddi [2001] T. Gehrmann and E. Remiddi, Numerical evaluation of harmonic polylogarithms, Comput. Phys. Commun. 141, 296 (2001), arXiv:hep-ph/0107173 .
  * Gehrmann and Remiddi [2002] T. Gehrmann and E. Remiddi, Numerical evaluation of two-dimensional harmonic polylogarithms, Comput. Phys. Commun. 144, 200 (2002), arXiv:hep-ph/0111255 [hep-ph] .
  * Vollinga and Weinzierl [2005] J. Vollinga and S. Weinzierl, Numerical evaluation of multiple polylogarithms, Comput. Phys. Commun. 167, 177 (2005), arXiv:hep-ph/0410259 .
  * Duhr and Dulat [2019] C. Duhr and F. Dulat, PolyLogTools — polylogs for the masses, JHEP 08, 135, arXiv:1904.07279 [hep-th] .
  * ’t Hooft and Veltman [1972] G. ’t Hooft and M. Veltman, Regularization and renormalization of gauge fields, Nucl. Phys. B44, 189 (1972).
  * Akyeampong and Delbourgo [1974] D. Akyeampong and R. Delbourgo, Anomalies via dimensional regularization, Nuovo Cim. A19, 219 (1974).
  * Larin [1993] S. Larin, The renormalization of the axial anomaly in dimensional regularization, Phys. Lett. B303, 113 (1993), arXiv:hep-ph/9302240 .
  * Lee [2012] R. N. Lee, Presenting LiteRed: a tool for the Loop InTEgrals REDuction, (2012), arXiv:1212.2685 [hep-ph] .
  * Argeri  _et al._ [2014] M. Argeri, S. Di Vita, P. Mastrolia, E. Mirabella, J. Schlenk, U. Schubert, and L. Tancredi, Magnus and Dyson series for master integrals, JHEP 03, 082, arXiv:1401.2979 [hep-ph] .