Hadronic Light by Light Corrections to the Muon Anomalous Magnetic Moment

Daniel Melo ${}^{a,\,}$111Corresponding author<EMAIL_ADDRESS>Edilson
Reyes<EMAIL_ADDRESS>and Raffaele Fazio
<EMAIL_ADDRESS>
a Departamento de Física, Universidad Nacional de Colombia, Bogotá D.C.,
Colombia.

b Departamento de Física, Universidad de Pamplona, Pamplona, Colombia.

###### Abstract

We review the Hadronic Light-by-Light (HLbL) contribution to the muon
anomalous magnetic moment. Upcoming measurements will reduce the experimental
uncertainty of this precision observable by a factor of four, thus breaking
the current balance with the theoretical prediction. A necessary step to
restore it is to decrease the HLbL contribution error, which implies a study
of the high-energy intermediate states that are neglected in dispersive
estimates. We focus on the maximally symmetric high-energy regime and in quark
loop approximation of perturbation theory we check the kinematic-
singularity/zero-free tensor decomposition of the HLbL amplitude.

## Introduction

The Standard Model (SM) is the current theoretical paradigm for particle
physics at its most fundamental level. This fact is rooted in the SM’s
mathematical consistency and specially in its highly accurate predictions for
precision experiments. In fact, one of the most precisely verified theoretical
predictions in the history of physics and the true triumph of quantum field
theory is the SM magnetic moment of the electron [1, 2, 3]
$\vec{\mu}=g\left(\frac{e}{2m}\right)\vec{S}$, being $m$ the electron mass and
$\vec{S}$ its spin operator. The so called anomalous part is expressed by the
quantity $a=\frac{g-2}{2}$, quantifying the deviation of the Land$\acute{e}$
factor from the classical value $g=2$, and is entirely due to quantum-
mechanical phenomena: the “cloud” of virtual particles with which the electron
is constantly interacting, slightly changes the way it interacts with a
classical magnetic field (see figure 1). Therefore, the measurement of the
anomalous part of a particle’s magnetic moment makes possible to test which
kind of other particles it interacts with and what is the strength of the
interaction. Consequently, this quantity is of the utmost interest for
theoretical physicists when testing the SM itself and also theories so called
Beyond the Standard Model (BSM). For the electron this anomalous part has been
computed to $O(\alpha^{5})$ in QED, for weak contributions the uncertainty is
$\sim 10^{-16}$ and for hadronic contributions it is $\sim 10^{-14}$ [3]. The
discrepancy with measurements is $8.76\cdot 10^{-13}$ [2] or $2.42$ times the
standard uncertainty (often represented as $\sigma$), which is still far
enough from the discovery threshold.

For the muon, the tension between the SM theoretical prediction for the
anomalous magnetic moment $a_{\mu}$ and its experimental measurement is
bigger. Therefore, it has attracted very much attention since the Brookhaven
National Laboratory (BNL) experiment results shed light on the issue in 2004
[4]. Currently, the tension stands at $4.2\sigma$ if the 2021 results from
Fermilab (FNAL) [5, 6] are taken into account in addition to the BNL ones:

$\displaystyle a_{\mu}^{BNL}$ $\displaystyle=116\;592\;089(63)\cdot
10^{-11}\;,$ (1) $\displaystyle a_{\mu}^{Fermilab}$
$\displaystyle=116\;592\;055(24)\cdot 10^{-11}\;,$ (2) $\displaystyle
a_{\mu}^{exp}$ $\displaystyle=116\;592\;059(22)\cdot 10^{-11}\;,$ (3)

As usual the combination of the two measurements is obtained from the
principle of maximum likelihood which, mathematically realized by the method
of least squares due to the Gaussian probability distribution, provides the
above weighted average.

In contrast, the most recent consensus SM prediction $a_{\mu}^{\text{SM}}$ is:

$a_{\mu}^{\text{SM}}=116\;591\;810(43)\cdot 10^{-11}$ (4)

which has been obtained by the “Muon $g-2$ Theory Initiative” and is described
in [7]. Consequently, the tension between the SM value and the measurement is:

$\displaystyle\Delta a_{\mu}$
$\displaystyle=a_{\mu}^{\text{exp}}-a_{\mu}^{\text{SM}}=249(48)\cdot
10^{-11}\text{ or }5.2\;\sigma\;.$ (5)

Although this discrepancy is already beyond the five sigma discovery
threshold, it is still not considered to be a sign of New Physics because
there are unresolved inconsistencies between the dispersive and lattice
estimations of the Hadronic Vacuum Polarization (HVP) contribution (see figure
2), which we will describe below.

From the three fundamental interactions considered in the SM, only the strong
force contribution currently has an uncertainty that is relevant with respect
to the tension’s value [7]. Therefore these strong contributions to $a_{\mu}$
are the main focus of the theoretical work towards reducing uncertainty.
Hadronic contributions affect $a_{\mu}$ in two ways: HVP and Hadronic Light by
Light Scattering (HLbL).

aMu  (5,5) i1,i2 o1,o2 phantomi2,v2,o2 fermioni1,v1,o1 photonv2,v1
arrow_len3mm (5,5) i1,i2 o1,o2 phantomi2,v2,o2 fermioni1,v1,o1 photonv2,v1
1cmv1 arrow_len3mm

Figure 1: Interaction of a fermion with a classical electromagnetic field at
tree level (left) vs. corrections due to virtual particles (right).

The HVP contribution is much larger than the HLbL one and moreover it can be
computed from experimental measurements in the well–known approach of
dispersive integrals [7]. More specifically, the contribution from HVP amounts
to $6845\times 10^{-11}$ and from HLBL it is $92\times 10^{-11}$. HVP is
essentially an hadronic correction to the photon propagator, then, because of
analyticity of Green functions and the unitarity of the theory, it can be
computed from the cross section of a virtual photon decaying into hadrons.
This cross section can be extracted from $e^{+}e^{-}\longrightarrow Hadrons$
data from detectors such as BaBar [8], KLOE [9], BESIII [10], SND [11] and
CMD–3 [12]. By the nature of the dispersive method, it is necessary to know
the hadron production cross section at different values of the center of mass
energy. This can be achieved either by directly changing the energy of the
$e^{-}$ and $e^{+}$ beams, called direct scan [13, 14], or by fixing it and
letting the (measured) initial–state radiation do the work of varying the
energy of the virtual photon which then decays into hadrons, called radiative
return [15, 16]. There are also alternative methods of measuring HVP by $\tau$
decay experiments [17] and by measuring the hadronic contribution to the
running of the fine structure constant $\alpha=e^{2}/4\pi$ from $\mu^{-}e^{-}$
elastic cross sections, called the MUonE project [18, 19, 20].

Before 2021, the main goals of the community for the HVP contribution were to
improve the uncertainty on dispersive estimates by solving tensions between
data sets, improving accountability of radiative corrections to the measured
cross sections and considering contributions from further channels. Even
though these goals still remain, the publication of the HVP contribution
estimate by the BMW lattice collaboration [21] upended the priorities. This
result, the first of its kind to have competitive uncertainties with respect
to the dispersive estimates, reduced the tension with experiment to 1.5
$\sigma$ when considered alone. Although this prediction is in tension with
dispersive estimates, it has been since confirmed by several other lattice
collaborations [22, 23, 24, 25, 26]. Interesting discussions regarding
comparison between dispersive and lattice estimates can be found in [27, 28,
29]. The most recent results from the CMD–3 experiment for $\pi\pi$ production
show a significant deviation from all other previous experimental results,
including CMD–2 [30]. Much like recent lattice HVP results, dispersive
estimates of HVP with $\pi\pi$ contributions coming only from this new set of
measurements also reduce the tension between the SM prediction and
$a_{\mu}^{\text{exp}}$ [29], which has further added to the confusion.
Consistency checks from analyticity and unitarity constraints on the pion
vector form factor have not shed much light on the discrepancies between
experiments [29].

HVP  (5,5) i1,i2 o1,o2 fermioni1,v4,v3,v1,o1 phantomi2,v2,o2 boson,photonv3,v2
photon,tension=0.01,label.side=leftv1,v5,v4 arrow_len3mm 0.9cmv5

Figure 2: Hadronic vacuum polarization contribution to the anomalous magnetic
moment of the muon. The blob contains only strongly interacting virtual
particles.

Now we go on to the HLbL scattering which gives a contribution (see figure 3)
of $92(18)\cdot 10^{-11}$. In contrast with HVP, the theoretical side of the
HLbL scattering computation had been much less understood until recently. The
added complexity is due to the fact that four currents are involved, instead
of only two. This introduces several difficulties. First of all, HLbL
scattering cannot be as cleanly related to $e^{+}e^{-}$ annihilation or other
experiments. Furthermore, the HLbL amplitude has a much more complex tensor
decomposition: it is a linear combination of 43 tensors, even after gauge
invariance constraints have been considered. Moreover, it is necessary to
expand this set to a redundant one with 54 elements in order to avoid
kinematic, meaning spurious or in general non-dynamical, singularities, that
spoil the dispersive approach. In the end, for the purpose of computing
$a_{\mu}^{\text{HLbL}}$, it is only necessary to know 7 of these scalar
coefficients, since the rest are related to them by crossing symmetry of
Mandelstam’s variables. Meanwhile, for HVP one initially has two tensor
structures which are then reduced to one due to gauge invariance. In fact,
this dispersion–fit tensor decomposition for HLbL was only recently found for
the first time [31, 32]. This multiplicity of scalar coefficients makes the
dispersive approach much more complex for HLbL than it is for HVP, because
each coefficient requires its own dispersive integral. In spite of this,
contribution of intermediate states including pseudoscalar poles, box
topologies and rescattering diagrams [7] have been successfully computed.
Particular applications with pions can be found in [33, 32, 34]. The
dispersive treatment of scalar and axial contributions is not yet in
satisfactory state, but significant progress has been made in that direction
lately [35, 36, 37, 38, 39, 40]. Before recent breakthroughs with the
dispersive method, the low energy regime of the HLbL scattering was studied
mostly with hadronic models, whose uncertainty was harder to assess. An
advantage on the computation of the HLbL contribution with respect to the HVP
one is that the first one appears at one further order of $\alpha$ than the
second one and thus its computation requires slightly less accuracy. Finally,
a common feature for both HLbL and HVP is that they are dominated by very
different degrees of freedom at low and high energies, namely, hadrons and
then quarks and gluons, respectively. The fact that HVP and HLbL amplitudes
enter the muon vertex as an insertion of one and two loops, respectively,
makes it necessary to properly “sew” the contributions from different
approaches at different kinematic regions.

HLbL (5,5) i o u fermioni,v1,v2,v3,o photon,tension=1v2,v4,u
photon,left=0.3,tension=1,label.side=left, label=$q_{1}$v1,v4
photon,right=0.3,tension=1,label.side=right, label=$q_{2}$v3,v4 0.9cmv4
$q_{4}\longrightarrow 0$u label=$q_{3}$,label.angle=135v2 arrow_len3mm

Figure 3: Hadronic light–by–light scattering contribution to the anomalous
magnetic moment of the muon. The blob contains only strongly interacting
virtual particles.

Much like with HVP, lattice estimations for HLbL have become competitive with
dispersive ones on recent years [41, 42, 43, 44, 45, 46]. Even though the
differences between the two are within uncertainty values, lattice estimates
also push the SM prediction towards the experimental value.

As mentioned previously, the HLbL amplitude at low photon virtualities (see
figure 3) was obtained from low–energy QCD models (scalar QED for the pion,
Nambu Jona-Lasinio model and vector meson dominance, for example) [47, 48] or,
more recently, from dispersive integrals on hadronic production from multiple
virtual photons [49]. Of course, these methods have a certain high energy
limit of validity, be it conceptual or practical. For the dispersive approach
it is the latter case. Extension to heavier intermediate states has been
hindered by a lack of data on the necessary subprocesses and the increasing
complexity of unitarity diagrams with multiple particles. Fortunately, heavy
intermediate states contributions are suppressed in dispersive integrals by a
narrower phase space and thus one can consider states up to certain mass and
still obtain a useful result. Nevertheless, to assess or reduce the
uncertainty coming from the neglected heavier states it is necessary to resort
to tools that complement, replace or evaluate the dispersive approach at high
energies. These tools are called short distance constraints (SDC). For
example, in the high energy regions of dispersive integrals, data for an
hadronic form factor can be replaced by the expression for its known
asymptotic behaviour. One can also evaluate the asymptotic behaviour of HLbL
scattering amplitude itself and use it to evaluate how well the set of
intermediate states considered in the dispersive approach resembles such
behaviour. For a finite number of intermediate states it is not possible to
completely mimic such behaviour [50, 51], but this fact can be used to measure
how well a set of intermediate states represent high energy contributions.
Such studies are a key complement of dispersive computations and play a
central role in uncertainty assessment [52, 53, 54]. There are two loop
momenta configurations that lead to a high energy regime in HLbL scattering:
$|q_{1}^{2}|\sim|q_{2}^{2}|\sim|q_{3}^{2}|\gg\Lambda_{QCD}^{2}$ and
$|q_{1}^{2}|\sim|q_{2}^{2}|\gg\Lambda_{QCD}^{2}$, where $q_{1}$, $q_{2}$ and
$q_{3}$ are the virtual photon momenta, $q_{4}$ is the real soft photon
momentum representing the electromagnetic external field (see figure 3) and
$\Lambda_{QCD}$ is the QCD hadronic regime threshold. The main purpose of this
research is to compute the HLbL scattering amplitude by the methods of
perturbative QCD in the regime in which the absolute values of the three
virtualities, $q_{1}^{2}$, $q_{2}^{2}$, $q_{3}^{2}$, are much larger than the
hadronic threshold. In particular, we perform an operator product expansion
(OPE) in an electromagnetic background for the HLbL scattering amplitude
following [55, 56, 57], but we fully harness the background field method to
provide an original and, in our view, more systematic framework, to include
the hadronic contributions in the same spirit of QCD sum rules. The main
result of the work is nevertheless the computation of the quark loop amplitude
that constitutes the leading contribution of the HLbL scattering amplitude at
high energies. Our computation can be considered an extension of the
literature’s result, because we obtain and study the full tensor structure of
the amplitude and obtain a complete series expansion of ligth quark mass
corrections up to arbitrary order. The computation is implemented using
original Mathematica scripts in combination with state–of–the–art packages,
FeynCalc [58, 59, 60] and MBConicHulls [61], for computations in high energy
physics.

##  1 $a_{\mu}$ in QFT

In this section we review the basics of the computation of the HLbL
contribution to $a_{\mu}$ with special focus on the dispersive approach and
the Mandelstam representation on which it is based. At the end of the section
we present a brief review of the low–energy contributions that are
well–understood so far.

###  1.1 Basics

The magnetic moment of a particle is defined through its scattering amplitude
on a classical magnetic field. More specifically, for a particle with spin $s$
and magnetic moment $\mu$ interacting with a classical magnetic field
$\boldsymbol{B}$, the matrix element of the interaction hamiltonian $H_{int}$
between an initial state $\psi_{\boldsymbol{p}\sigma}$ with momentum
$\boldsymbol{p}$ and spin projection $\sigma$ and final state
$\psi_{\boldsymbol{p}^{\prime}\sigma^{\prime}}$ is:

$\displaystyle\langle\psi_{\boldsymbol{p}^{\prime}\sigma^{\prime}}|H_{int}|\psi_{\boldsymbol{p}\sigma}\rangle$
$\displaystyle=-\frac{\mu}{s}(\boldsymbol{J}^{(s)})_{\sigma^{\prime}\sigma}\cdot\boldsymbol{B}\;\delta^{3}(\boldsymbol{p^{\prime}}-\boldsymbol{p})\times
2m\;,$ (6)

where $\delta^{3}$ represents the Dirac delta in three dimensions, $m$ stands
for the particle’s mass and $\boldsymbol{J}^{(s)}$ is the little group
generator associated to a massive particle of spin $s$. The factor $2m$
appears only due to the relativistic normalization of the states.

For a relativistic charged particle, the corresponding matrix element is:

$\displaystyle\langle\psi_{\boldsymbol{p}^{\prime}\sigma^{\prime}}|H_{int}|\psi_{\boldsymbol{p}\sigma}\rangle$
$\displaystyle=-e_{q}j_{\mu}A^{\mu}\;,$ (7)

where $A_{\mu}$ is the classical electromagnetic potential, $j_{\mu}$ is the
matrix element of the particle’s current operator and $e_{q}$ represents its
electric charge. For the muon we have $s=1/2$, $e_{q}=-e$444$e$ represents the
absolute value of the electric charge of the electron. and:

$\displaystyle j_{\mu}(x)$
$\displaystyle=e^{i(p-p^{\prime})x}\langle\mu^{-}_{\boldsymbol{p}^{\prime}\sigma^{\prime}}|J^{\mu}(0)|\mu^{-}_{\boldsymbol{p}\sigma}\rangle=e^{i(p-p^{\prime})x}\overline{u}_{\boldsymbol{p}^{\prime}\sigma^{\prime}}\Gamma_{\mu}(p^{\prime},p)u_{\boldsymbol{p}\sigma}\;,$
(8)

where $J^{\mu}$ represents the electromagnetic current Heisenberg operator of
the muon, $\Gamma_{\mu}$ (in QFT) is the amplitude of the full on–shell vertex
diagram to the right of Figure 1 and $u_{\boldsymbol{p}s}$ and
$\overline{u}_{\boldsymbol{p}^{\prime}s^{\prime}}$ are the spinors associated
to the incoming and outgoing muon, respectively.

Considering the fact that $j_{\mu}$ has to behave as a Lorentz vector and it
may contain Dirac matrices, then $\Gamma_{\mu}$ must be a linear combination
of the four–momenta, the Levi–Civita symbol $\epsilon^{\mu\nu\lambda\rho}$ and
Dirac bilinears, which are listed in the first column of table 1. The complete
set of independent structures that can be built including index contractions
is written in the second column table 1, where the total momentum
$P^{\mu}\equiv p^{\mu}+p^{{}^{\prime}\mu}$ and the transferred momentum
$q^{\mu}\equiv p^{{}^{\prime}\mu}-p^{\mu}$ have been used. However, such set
can be greatly reduced using the on–shell character of the spinors in (8)
(i.e. Dirac’ equation) and gauge invariance.

Dirac matrices’ basis element | Available structures
---|---
1 |  $P^{\mu}$, $q^{\mu}$
$\gamma^{\mu}$ |  $\gamma^{\mu}$, $q^{\mu}\not{q}$, $q^{\mu}\not{P}$, $P^{\mu}\not{q}$, $P^{\mu}\not{P}$, $\epsilon^{\mu\nu\lambda\rho}\gamma_{\nu}P_{\lambda}q_{\rho}$
$\gamma^{5}$ |  $P^{\mu}\gamma^{5}$, $q^{\mu}\gamma^{5}$
$\gamma^{\mu}\gamma^{5}$ |  $\gamma^{\mu}\gamma^{5}$, $q^{\mu}\not{q}\gamma^{5}$, $q^{\mu}\not{P}\gamma^{5}$, $P^{\mu}\not{q}\gamma^{5}$, $P^{\mu}\not{P}\gamma^{5}$, $\epsilon^{\mu\nu\lambda\rho}\gamma_{\nu}\gamma^{5}P_{\lambda}q_{\rho}$
$\sigma_{\mu\nu}$ |  $P^{\alpha}q^{\beta}\sigma_{\alpha\beta}P^{\mu}$, $P^{\alpha}q^{\beta}\sigma_{\alpha\beta}q^{\mu}$, $\sigma^{\mu\nu}P_{\nu}$, $\sigma^{\mu\nu}q_{\nu}$, $\epsilon^{\mu\nu\lambda\rho}\sigma_{\nu\lambda}P_{\rho}$, $\epsilon^{\mu\nu\lambda\rho}\sigma_{\nu\lambda}q_{\rho}$
Table 1: A priori available structures in the covariant decomposition of the
muon electromagnetic on–shell vertex $\Gamma^{\mu}$ (see equation (8) and
figure 1). We use $\sigma_{\mu\nu}=\frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$

Indeed the most general tensor structure of $\Gamma_{\mu}(p^{\prime},p)$ for
on–shell spinors is:

$\Gamma_{\mu}(p^{\prime},p)=A_{1}(q^{2})\gamma_{\mu}+P^{\mu}A_{2}(q^{2})+\left(\gamma^{\mu}-\frac{2mq^{\mu}}{q^{2}}\right)\gamma^{5}A_{3}(q^{2})+P^{\mu}\gamma^{5}A_{4}(q^{2})\;.$
(9)

By using the Gordon’s identity:

$\displaystyle\overline{u}_{\boldsymbol{p}^{\prime}\sigma^{\prime}}\gamma^{\mu}u_{\boldsymbol{p}\sigma}$
$\displaystyle=\frac{1}{2m}\overline{u}_{\boldsymbol{p}^{\prime}\sigma^{\prime}}\\{P^{\mu}+i\sigma^{\mu\nu}q_{\nu}\\}u_{\boldsymbol{p}\sigma}\;,$
(10) $\displaystyle\Gamma_{\mu}(p^{\prime},p)$
$\displaystyle=F_{1}(q^{2})\gamma_{\mu}+i\sigma^{\mu\nu}\frac{q_{\nu}}{2m}F_{2}(q^{2})+(\gamma^{\mu}-\frac{2mq^{\mu}}{q^{2}})\gamma^{5}F_{3}(q^{2})+\sigma^{\mu\nu}\frac{q_{\nu}}{2m}\gamma^{5}F_{4}(q^{2})\;.$
(11)

This convention is helpful because
$\sigma_{\mu\nu}\equiv\frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$ is the generator
of Lorentz transformations for covariant wave functions of Dirac fermions and
therefore $\sigma_{ij}$ generates rotations and little group transformations.
Thus the parallel with (6) becomes straightforward. In (11) $F_{1}$, $F_{2}$,
$F_{3}$ and $F_{4}$ are Lorentz invariant coefficients, also called “form
factors”. The first two are associated to parity conserving contributions and
are also known as electric and magnetic form factor, respectively. On the
other hand, $F_{3}$ and $F_{4}$ are related to parity violating and CP
violating contributions, respectively, and are also known as anapole moment
and electric dipole moment. For a pure on-shell derivation of the structure of
the electromagnetic vertex including, after an appropriate analytic
continuation of momenta, also the photon to be on-shell we refer to [62],
[63]. Since the anomalous magnetic moment is related to a non–relativistic
interaction, it is necessary to evaluate the muon vertex in the limit of zero
exchanged momentum, that is $q\rightarrow 0$. In such limit we have
$F_{1}(0)=1$ in order to define $e$ as the physical electric charge measured
in the interaction with a classical Coulomb field. We also have $F_{3}(0)=0$
in this limit. On the other hand, $F_{2}(0)$ and $F_{4}(0)$ are not
constrained. Thus, in the limit of zero exchanged momentum and a slowly
varying magnetic field (7) and (8) become:

$\displaystyle\langle\psi_{\boldsymbol{p}^{\prime}\sigma^{\prime}}|H_{int}|\psi_{\boldsymbol{p}\sigma}\rangle$
$\displaystyle=2m\times\frac{e\,}{m}(1+F_{2}(0))(\boldsymbol{J}^{(1/2)})_{\sigma^{\prime}\sigma}\cdot\boldsymbol{B}\;\delta^{3}(\boldsymbol{p^{\prime}}-\boldsymbol{p})\;,$
(12)

which implies:

$\displaystyle\mu$ $\displaystyle=\frac{e}{2m}(1+F_{2}(0))\;.$ (13)

At tree level we have $\Gamma_{\mu}=\gamma_{\mu}$ and thus
$\mu=\frac{e_{q}}{2m}$, which agrees with Pauli’s equation and Dirac’s
equation in the non–relativistic limit. Then, quantum corrections to this
classical value can be singled out by the gyromagnetic factor $g$:

$\displaystyle\mu\equiv g\frac{e_{q}}{2m}s\implies
a\equiv\frac{1}{2}(g-2)=F_{2}(0)\;,$ (14)

where $a$ is called the anomalous part of the magnetic moment of such
particle. For the muon we will use the symbols $a_{\mu}$ and $g_{\mu}$.

Once it is clear which part of the on–shell vertex in (8) actually contributes
to $\mu$, it is convenient to project it out. The projector needed for the
magnetic form factor is [64]:

$\displaystyle P_{2}^{\mu}$
$\displaystyle\equiv-\frac{m^{2}}{q^{2}(q^{2}-4m^{2})}\Big{(}\gamma^{\mu}+\frac{q^{2}+2m^{2}}{m(q^{2}-4m^{2})}P^{\mu}\Big{)}\;,$
(15)

which is to be used in the following way:

$\displaystyle F_{2}(q^{2})$
$\displaystyle=\operatorname{Tr}\\{(\not{p}+m)P_{2\mu}(\not{p}^{\prime}+m)\Gamma^{\mu}\\}\;,$
(16)

Note that there is a divergent $1/q^{2}$ factor in $P_{2\mu}$. Although
expanding $\Gamma^{\mu}$ around $q=0$ is inevitable, we can truncate the
expansion at first order:

$\displaystyle\Gamma^{\mu}(q^{2})$
$\displaystyle=\Gamma^{\mu}(0)+q_{\nu}\underbrace{\partial^{\nu}\Gamma^{\mu}(q^{2})|_{q=0}}_{\equiv\Gamma^{\mu\nu}}\;.$
(17)

To prove that, we must first reorganize $P_{2}$ using a slightly different
version of the Gordon identity:

$\displaystyle(\not{p}^{\prime}+m)\gamma^{\mu}(\not{p}+m)$
$\displaystyle=\frac{1}{2m}(\not{p}^{\prime}+m)\\{P^{\mu}+i\sigma^{\mu\nu}q_{\nu}\\}(\not{p}+m)\;,$
(18) $\displaystyle(\not{p}+m)P_{2}^{\mu}(\not{p}^{\prime}+m)$
$\displaystyle=(\not{p}+m)\frac{-m}{2(q^{2}-4m^{2})}\Big{(}-i\sigma^{\mu\nu}\frac{q_{\nu}}{q^{2}}+\frac{3}{(q^{2}-4m^{2})}P^{\mu}\Big{)}(\not{p}^{\prime}+m)\;.$
(19)

We see that the divergent term in the projector is of order $1/q$, so the
conclusion follows.

By inserting (17) as well as (19) into (16) with the substitutions
$p^{\prime}=\frac{1}{2}(P+q)$ and $p=\frac{1}{2}(P-q)$ we evaluate the
corresponding expression at $q=0$ obtaining after some algebra:

$\begin{split}a_{\mu}&=\lim_{q\rightarrow
0}\frac{im}{2(q^{2}-4m^{2})}\times\\\
&\operatorname{Tr}\Big{\\{}\Big{(}-\frac{1}{2}\not{q}\sigma_{\mu\nu}\frac{q^{\nu}}{q^{2}}\Big{[}\frac{\not{P}}{2}+m\Big{]}+\Big{[}\frac{\not{P}}{2}+m\Big{]}\sigma_{\mu\nu}\frac{q^{\nu}}{q^{2}}\Big{[}\frac{\not{P}}{2}+m\Big{]}+\frac{1}{2}\Big{[}\frac{\not{P}}{2}+m\Big{]}\sigma_{\mu\nu}\frac{q^{\nu}}{q^{2}}\not{q}\\}\Big{)}\Gamma^{\mu}(0)\Big{\\}}\\\
&+\lim_{q\rightarrow
0}\frac{im}{2(q^{2}-4m^{2})}\operatorname{Tr}\Big{\\{}\Big{[}\frac{\not{P}}{2}+m\Big{]}\sigma_{\mu\nu}\frac{q^{\nu}q_{\beta}}{q^{2}}\Big{[}\frac{\not{P}}{2}+m\Big{]}\Gamma^{\mu\beta}\Big{\\}}-\frac{3}{8m^{2}}\operatorname{Tr}\Big{\\{}p_{\mu}(m+\not{p})\Gamma^{\mu}(0)\Big{\\}}\;,\end{split}$
(20)

where we have explicitly evaluated the limit when possible. We still have
divergent terms together with tensor dependence on $q^{\mu}$. To get rid of
the latter we will take advantage of the scalar character of $F_{2}(q^{2})$ to
perform Lorentz transformations on $q^{\mu}$ before taking the $q\rightarrow
0$ limit.

In particular, we can perform spatial rotations, thus, we carry out an angular
average over the spatial components of $q^{\mu}$ taking $P$ as reference. The
results are:

$\int\frac{d\Omega}{4\pi}q^{\mu}=0\;,\hskip
28.45274pt\int\frac{d\Omega}{4\pi}q^{\mu}q^{\nu}=\frac{q^{2}}{3}\Big{(}g^{\mu\nu}-\frac{P^{\mu}P^{\nu}}{P^{2}}\Big{)}\;.$
(21)

The first result is obvious by the oddness of the integrand. The tensor
structure of the second result is evident from Lorentz covariance and the
scalar factors can be obtained straightforwardly by computing the Lorentz
trace on both sides.

After inserting the angular averages (21) inside (20) we obtain:

$a_{\mu}=\operatorname{Tr}\Big{\\{}\Big{(}\frac{1}{12}\gamma_{\mu}-\frac{1}{4}\frac{p_{\mu}}{m}-\frac{1}{3}\frac{1}{m^{2}}p_{\mu}\not{p}\Big{)}\Gamma^{\mu}(0)+\frac{1}{48m}(\not{p}+m)[\gamma_{\mu},\gamma_{\beta}](\not{p}+m)\Gamma^{\mu\beta}\Big{\\}}\;,$
(22)

which is the direct relation between $a_{\mu}$ and Feynman amplitudes
($\Gamma^{\mu}$) that we were looking for and opens the path to specialize the
result to particular topologies of diagrams and, in particular, to the HLbL
one.

###  1.2 Specializing $a_{\mu}$ to HLbL scattering amplitudes

In this section we specialize the result obtained in previous one to compute
$a_{\mu}$ from HLbL scattering amplitudes, whose diagrams are shown in Figure
3. The term “light–by–light” makes reference to the subdiagram appearing in
figure 3, which has four external photons (three virtual and attached to the
muon line and one representing and external field). The term “hadronic” is due
to the fact that only strongly interacting particles (quarks and gluons) or
hadrons (mesons and various resonances) are allowed to appear in the blob of
figure 3, either as virtual exchanged particles or as poles of the amplitude,
for the HLbL contributions.

The first step to specialize $a_{\mu}$ to HLbL is to isolate the appropriate
subdiagram amplitudes from the muon electromagnetic vertex ones. Making use of
the Feynman rules for QED it is possible to read the result off the Feynman
diagram in figure 3:

$\displaystyle-e\overline{u}_{\boldsymbol{p}^{\prime}\sigma^{\prime}}\Gamma^{\mu_{4}}u_{\boldsymbol{p}\sigma}$
$\displaystyle=\overline{u}_{\boldsymbol{p}^{\prime}\sigma^{\prime}}\int\frac{d^{4}q_{1}}{(2\pi)^{4}}\int\frac{d^{4}q_{2}}{(2\pi)^{4}}\frac{-i}{q_{1}^{2}}\frac{-i}{q_{2}^{2}}\frac{-i}{q_{3}^{2}}(-ie\gamma_{\mu_{1}})i\frac{\not{p}^{\prime}+\not{q}_{1}+m}{(p^{\prime}+q_{1})^{2}-m^{2}}(-ie\gamma_{\mu_{3}})$
$\displaystyle\hskip 28.45274pt\times
i\frac{\not{p}-\not{q}_{2}+m}{(p-q_{2})^{2}-m^{2}}(-ie\gamma_{\mu_{2}})\times
e^{4}\,\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(q_{1},q_{2},q_{3})u_{\boldsymbol{p}\sigma}\;,$
(23)

where $e^{4}\,\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ represents the amplitude of
the hadronic blob inside figure 3:

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ $\displaystyle=-i\int
d^{4}x\int d^{4}y\int
d^{4}ze^{-i(q_{1}x+q_{2}y+q_{3}z)}\langle\Omega|J^{\mu_{1}}_{s}(x)J^{\mu_{2}}_{s}(y)J^{\mu_{3}}_{s}(z)J^{\mu_{4}}_{s}(0)|\Omega\rangle\;.$
(24)

$|\Omega\rangle$ represents the QCD vacuum and $J_{s}$ stands for the
electromagnetic current of strongly interacting particles. In the literature,
$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ is referred to as “fourth rank vacuum
polarization tensor” [65, 66, 67] or “HLbL tensor” [49, 55]. Introducing this
new convention into (23) we obtain:

$\Gamma^{\mu_{4}}_{\text{HLbL}}=-e^{6}\int\frac{d^{4}q_{1}}{(2\pi)^{4}}\int\frac{d^{4}q_{2}}{(2\pi)^{4}}\frac{\gamma_{\mu_{1}}}{q_{1}^{2}}\frac{\not{p}^{\prime}+\not{q}_{1}+m}{(p^{\prime}+q_{1})^{2}-m^{2}}\frac{\gamma_{\mu_{3}}}{q_{3}^{2}}\frac{\not{p}-\not{q}_{2}+m}{(p-q_{2})^{2}-m^{2}}\frac{\gamma_{\mu_{2}}}{q_{2}^{2}}\,\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}\;.$
(25)

The next step towards finding the HLbL contribution to $a_{\mu}$, labelled
$a_{\mu}^{\text{HLbL}}$, is to insert (25) into (22).

The $\Gamma^{\mu}(0)$ term vanishes for HLbL, as can be deduced from the
analysis presented in [68], which concludes that the cross section of a
process in the limit in which an external photon becomes soft is equal to a
sum of terms proportional to the amplitude of the process without the soft
photon and its derivative plus vanishing contributions proportional to the
soft photon momentum. In the context of HLbL scattering, the previous
statement reads:

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}(q_{1},q_{2},q_{3})$
$\displaystyle\sim\Pi^{\mu_{1}\mu_{2}\mu_{3}}A^{\mu_{4}}+\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}}B^{\mu_{4}}+O(q_{4})\;,$
(26)

where $\Pi^{\mu_{1}\mu_{2}\mu_{3}}$ represents the three–photon scattering
amplitude, $A_{\mu_{4}}$ and $B_{\mu_{4}}$ are two vectors of order
$O(1/q_{4})$ and $O(q_{4}^{0})$, respectively, and $\partial$ represents
derivative with respect to some kinematic variable of the problem.
$\Pi^{\mu_{1}\mu_{2}\mu_{3}}$ vanishes due to Furry’s theorem and therefore
$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ vanishes (at least) linearly in the
static field limit555It is not possible to arrive at this conclusion using
Weinberg’s eikonal factor since the external soft photon may be emitted by an
internal line in the strongly–interacting blob..

With respect to $\Gamma^{\mu\alpha}_{\text{HLbL}}$, from (25) we find:

$\Gamma^{\mu_{4}\nu_{4}}=e^{6}\int\frac{d^{4}q_{1}}{(2\pi)^{4}}\int\frac{d^{4}q_{2}}{(2\pi)^{4}}\frac{\gamma_{\mu_{1}}}{q_{1}^{2}}\frac{\not{p}+\not{q}_{1}+m}{(p^{\prime}+q_{1})^{2}-m^{2}}\frac{\gamma_{\mu_{3}}}{q_{3}^{2}}\frac{\not{p}-\not{q}_{2}+m}{(p-q_{2})^{2}-m^{2}}\frac{\gamma_{\mu_{2}}}{q_{2}^{2}}\,\partial^{\mu_{4}}\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}\Big{|}_{q_{4}\rightarrow
0}\;,$ (27)

where we have used the antisymmetry between $\mu_{4}$ and $\nu_{4}$ of
$\partial^{\mu_{4}}\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}|_{q_{4}\rightarrow 0}$,
which can be deduced by differentiating the Ward identity. Terms proportional
to $q_{4}^{\nu_{4}}\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ were disregarded in the
static limit. Finally, turning back to $a_{\mu}^{\text{HLbL}}$ one obtains:

$\displaystyle a_{\mu}^{\text{HLbL}}$
$\displaystyle=\frac{e^{6}}{48m}\int\frac{d^{4}q_{1}}{(2\pi)^{4}}\int\frac{d^{4}q_{2}}{(2\pi)^{4}}\frac{1}{q_{1}^{2}}\frac{1}{q_{2}^{2}}\frac{1}{q_{3}^{2}}\frac{1}{(p+q_{1})^{2}-m^{2}}\frac{1}{(p-q_{2})^{2}-m^{2}}\,\frac{\partial}{\partial
q_{4\mu_{4}}}\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}\Big{|}_{q_{4}\rightarrow 0}$
$\displaystyle\hskip
28.45274pt\times\operatorname{Tr}\Big{\\{}(\not{p}+m)[\gamma_{\mu_{4}},\gamma_{\nu_{4}}](\not{p}+m)\gamma_{\mu_{1}}(\not{p}+\not{q}_{1}+m)\gamma_{\mu_{3}}(\not{p}-\not{q}_{2}+m)\gamma_{\mu_{2}}\Big{\\}}\;.$
(28)

There are only three steps left to obtain $a_{\mu}^{\text{HLbL}}$: (i) compute
the Dirac trace, (ii) compute
$\partial^{\mu_{4}}\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}|_{q_{4}\rightarrow 0}$
and (iii) compute the two–loop integral. The trace can be performed
straightforwardly. On the other hand, the computation of the HLbL amplitude is
very complex and it is therefore necessary to study it in depth before
advancing further.

###  1.3 Dispersive computation of the HLbL amplitude

In the previous section we expressed $a_{\mu}$ in terms of the HLbL scattering
amplitude $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ and in this section we will
present the dispersive approach for computing it.

Since $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ appears inside a two loop integral
on $q_{1}$ and $q_{2}$, it is necessary to compute it at different energy
regions involving perturbative and non–perturbative regimes. This is due to
asymptotic freedom, which invalidates perturbation theory at energy scales
below $\Lambda_{QCD}\sim$1 GeV because the coupling $\alpha_{s}$ approaches
$1$. Non–perturbative contributions give the bulk of $a_{\mu}^{\text{HLbL}}$,
but high energy studies are important for error estimation, so it is
convenient to perform computations in both regimes in a unified framework.

QFT in the lattice and dispersive integrals are two of very few tools that
allow for the computation of amplitudes in non–perturbative regimes. The first
one tries to solve the QFT equations in a finite spacetime cube of side length
$L$ with discrete Euclidean spacetime coordinates of spacing $a$. The
observables of interest are then computed for different values of large $L$
and $1/a$ and these results are then extrapolated to $L,1/a\rightarrow\infty$
in order to recover the standard QCD results. Using a very different
perspective, the dispersive approach [49, 31, 32, 33] relies on the
analyticity of scattering amplitudes and the unitarity of the S–matrix
(probability conservation) to relate the amplitudes of a process with the
cross sections of its sub–processes, which are fitted to data. Lattice
computations have high numerical complexity due to the very large number of
degrees of freedom that appear when both the size $L$ of the system and its
resolution $1/a$ become large, as needed to reduce the uncertainty produced by
the extrapolation to $L,1/a\rightarrow\infty$. Therefore, even for the simpler
case of HVP only till very recently have its results become competitive with
the dispersive ones in terms of uncertainty [69, 21, 45]. Furthermore, for
HLbL the first lattice computations are relatively new and they still are not
competitive with dispersive ones. Although the dispersive approach for the
computation of the HLbL contribution to $a_{\mu}$ also has its drawbacks, the
main one relating to the Mandelstam representation of
$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ have been recently overcome. This has
allowed to obtain the most reliable accounts of $a_{\mu}^{\text{HLbL}}$ in
recent years [7]666Although the approach proposed in [49, 31, 32] has been the
standard scheme in recent estimations of $a_{\mu}^{\text{HLbL}}$, there exist
alternate dispersive frameworks. For example, in [70, 71, 72, 73, 74] a
dispersive equation is applied directly to the magnetic form factor $F_{2}$
instead of the HLbL Feynman amplitude..

As mentioned at the beginning of this subsection, in this work we focus on the
dispersive approach. It is based on four fundamental pillars: unitarity of the
$S$–matrix, the Sugawara–Kanazawa theorem for functions of a complex variable
and the Schwarz reflection theorem. We will review such pillars in order.

####  1.3.1 Unitarity of the S–matrix

A key concept in relativistic quantum theories of fundamental
interactions777Not any quantum theory conserves probability, non–relativistic
systems with absorptive potentials are an example. is conservation of
probability. For transition rates, it is equivalent to the unitarity of the
S–matrix. Let us explore the consequences of such feature for the transition
matrix $T$:

$\displaystyle S^{\dagger}$
$\displaystyle=S^{-1}\implies(1+iT)(1-iT^{{\dagger}})=1\implies
TT^{{\dagger}}=i(T^{{\dagger}}-T)\;.$ (29)

If we evaluate a certain matrix element of $S$ and insert a complete set of
momentum eigenstates in the last equation we obtain [75]:

$\displaystyle 2\text{Im}\mathcal{M}(i\rightarrow f)$
$\displaystyle=\displaystyle\sum_{n}\Big{(}\Pi^{n}_{i=1}\int\frac{d^{3}\boldsymbol{q}_{i}}{(2\pi)^{3}}\frac{1}{2E_{i}}\Big{)}\mathcal{M}^{*}(f\rightarrow\\{q_{i}\\})\mathcal{M}(i\rightarrow\\{q_{i}\\})\,(2\pi)^{4}\delta^{4}(P_{i}-\sum_{i}q_{i})\;,$
(30)

where $\mathcal{M}$ stands for a Feynman amplitude, $i$ represents the initial
state, $f$ the final one and ${q_{i}}$ the on–shell intermediate–states, which
come from the insertion of the identity resolved in terms of the momentum
eigenstates. If the initial and final states are the same, like for HVP (see
figure 2), then we obtain the optical theorem.

At this point (30) does not seem to be of much help. First, it can only be
applied in principle to the on–shell electromagnetic vertex amplitude.
Secondly, we only have the imaginary part of
$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ and we need it in full to compute
observables such as $a_{\mu}^{\text{HLbL}}$ in (22). Additionally, we are
tasked to compute one amplitude in terms of infinitely many and infinitely
complex different subamplitudes. The first issue can be fixed by simplifying
terms in both sides of the equation such that it applies just as well to the
HLbL subdiagrams. To deal with the second one we need to make use of the
Sugawara–Kanazawa theorem, which reconstructs a complex variable function
based on its poles and branch cuts. It will also show us why we can consider
only intermediate states up to certain energy scale in (30) and still get
meaningful predictions.

####  1.3.2 Sugawara–Kanazawa theorem and Schwarz reflection identity

Consider a function of a complex variable $z$ with (possibly) two branch cuts
along the real axis: one to the right starting at $c_{1}$ and extending
(possibly) to positive infinity and one to the left starting at $-c_{2}$ and
extending (possibly) to negative infinity. Based on the following three
requirements:

  * •

$f(z)$ has finite limits in the positive real infinity direction above and
below the right–hand cut.

  * •

The limit of $f(z)$ in the negative real infinity direction above and below
the left–hand cut exists.

  * •

If $f(z)$ is divergent in certain infinite direction, such a divergence is
weaker than a polynomial with finite power $N$ such that $N\geq 1$.

then the Sugawara–Kanazawa theorem [76] claims that $f(z)$ may be represented
as888There is also another result to this theorem that essentially claims that
$f(z)$ has the same limit at infinity in any direction with positive
(negative) imaginary part as it has along and above (under) the right (left)
cut.:

$\displaystyle
f(z)=\displaystyle\sum_{i}\frac{R_{i}}{z-x_{i}}+\frac{1}{\pi}\Big{(}\int_{c_{1}}^{\infty}+\int_{-\infty}^{-c_{2}}\Big{)}\frac{\Delta_{x}f(x)}{x-z}\;dx+\lim_{x\rightarrow\infty}\overline{f}(x)\;,$
(31)
$\displaystyle\Delta_{x}f(x)=\frac{1}{2i}\\{f(x+i\epsilon)-f(x-i\epsilon)\\}\;,\hskip
28.45274pt\overline{f}(x)=\frac{1}{2}\\{f(x+i\epsilon)+f(x-i\epsilon)\\}\;,$

where $R_{i}$ represents the residue of $f(z)$ in $x_{i}$ which lies on the
real interval $[-c_{2},c_{1}]$. The two integrals in (31) are performed along
the real axis. This representation of $f(z)$ is usually called “dispersion
relation”. The last term is referred to as the “subtraction constant” and
accounts for possible divergences of $f(z)$ in the infinity, which enter the
equation as the contribution of the circumference of a Cauchy integration path
at infinity.

The HLbL process may be regarded as a function of two of the usual Mandelstam
variables $s$, $t$ and $u$ for two–two scattering. Therefore, the
Sugawara–Kanazawa theorem has to be applied twice; one for each independent
variable, leaving the rest constant. However, the procedure is to be applied
to a function that unifies all three channels’ amplitudes. The resulting
double dispersive integral is known as Mandelstam representation [77]. The
first step to build it is to write a single dispersive representation for,
say, $s$, which we will consider to be unsubtracted for simplicity. Since we
are expecting $\mathcal{M}$ to “contain” the amplitudes for the three
channels, it must be invariant under crossing. As such, given that $t$ has a
fixed value, we expect to have:

$\displaystyle\mathcal{M}(s,t)$
$\displaystyle=\displaystyle\sum_{i}\frac{R_{i}^{s}(t)}{s-x^{s}_{i}}+\frac{1}{\pi}\int_{c_{1}}^{\infty}\frac{\Delta_{s^{\prime}}\mathcal{M}(s^{\prime},t)}{s^{\prime}-s}\;ds^{\prime}+\frac{1}{\pi}\int_{c_{1}}^{\infty}\frac{\Delta_{u^{\prime}}\mathcal{M}(u^{\prime},t)}{u^{\prime}-u}\;du^{\prime}\;.$
(32)

Now we perform an analytic continuation on $t$ to whichever value we require.
Via interaction form factors the $t$ dependence of the residues is usually
well–known, but not for $\Delta\mathcal{M}$. To deal with this we apply (31)
once again, but this time for $\Delta_{s}\mathcal{M}$ with fixed $s^{\prime}$:

$\Delta_{s^{\prime}}\mathcal{M}(s^{\prime},t)=\frac{1}{\pi}\int_{c_{3}(s^{\prime})}^{\infty}\frac{\Delta_{t^{\prime}}\Delta_{s^{\prime}}\mathcal{M}(s^{\prime},t^{\prime})}{t^{\prime}-t}\;dt^{\prime}+\frac{1}{\pi}\int_{c_{3}(s^{\prime})}^{\infty}\frac{\Delta_{u^{\prime}}\Delta_{s^{\prime}}\mathcal{M}(s^{\prime},u^{\prime})}{u^{\prime}-\overline{u}}\;du^{\prime}$
(33)

where $\overline{s}$ is the Mandelstam variable associated to $u^{\prime}$ and
$t$, while $\overline{u}$ is the Mandelstam variable associated to
$s^{\prime}$ and $t$. An similar equation applies for $\Delta_{u}\mathcal{M}$.
Inserting this into the single dispersion relation we obtain:

$\displaystyle\mathcal{M}(s,t)$
$\displaystyle=\frac{1}{\pi^{2}}\int_{c_{1}}^{\infty}ds^{\prime}\int_{c_{3}(s^{\prime})}^{\infty}dt^{\prime}\frac{\Delta_{t^{\prime}}\Delta_{s^{\prime}}\mathcal{M}(s^{\prime},t^{\prime})}{(s^{\prime}-s)(t^{\prime}-t)}+\frac{1}{\pi^{2}}\int_{c_{1}}^{\infty}ds^{\prime}\int_{c_{3}(s^{\prime})}^{\infty}du^{\prime}\frac{\Delta_{u^{\prime}}\Delta_{s^{\prime}}\mathcal{M}(s^{\prime},u^{\prime})}{(s^{\prime}-s)(u^{\prime}-u)}$
$\displaystyle+\frac{1}{\pi^{2}}\int_{c_{1}}^{\infty}du^{\prime}\int_{c_{3}(u^{\prime})}^{\infty}dt^{\prime}\frac{\Delta_{t^{\prime}}\Delta_{s^{\prime}}\mathcal{M}(t^{\prime},u^{\prime})}{(u^{\prime}-u)(t^{\prime}-t)}+\displaystyle\sum_{i}\frac{R_{i}^{s}(t)}{s-x^{s}_{i}}\;.$
(34)

Some terms were simplified by noting that
$s^{\prime}-\overline{s}=s^{\prime}-s+u^{\prime}-u=u^{\prime}-\overline{u}$.

The $\Delta_{i}\Delta_{j}\mathcal{M}$ are called double spectral functions and
can be obtained from the optical theorem with help from the Schwarz reflection
principle. It states that if a function $f(z)$ of a complex variable $z$ is
real along certain finite segment $\Gamma$ of the real axis, then
$f^{*}(z)=f(z^{*})$ in a domain $D$ of the $z$ complex plane that contains
$\Gamma$ and in which $f(z)$ is analytic. This is always the case for any
amplitude because at a sufficiently low center–of–mass (CM) energy no
intermediate state is allowed and thus the amplitude becomes real.

On the other hand, this theorem implies the existence of a discontinuity
across the physical region of the real axis for each kinematic variable
because there is always at least one possible intermediate state for any
process: the initial one. Let $z$ be a kinematic variable of an amplitude
$\mathcal{M}$ in a physical region, then:

$2i\text{Im}\mathcal{M}(z+i\epsilon)=\mathcal{M}(z+i\epsilon)-\mathcal{M}^{*}(z+i\epsilon)=\mathcal{M}(z+i\epsilon)-\mathcal{M}(z-i\epsilon)=\Delta\mathcal{M}(z)\;.$
(35)

This result relates the Mandelstam representation with (30). We can see that
the constant $c_{1}$ of the dispersive integral is in fact the CM–frame energy
of the lightest multiparticle intermediate state. One–particle intermediate
states correspond to poles. Analogously, this helps us understand the meaning
of $c_{3}$. For a physical value of $s$ and $t$, $\Delta\mathcal{M}$ is just
$\text{Im}\mathcal{M}$ and is hence real. However, if we analytically continue
$\Delta\mathcal{M}$ beyond the physically permissible boundaries of $t$, it
may (and does) become complex, betraying the existence of a discontinuity.
$c_{2}(s)$ is the point where that happens. Note that this is true even if $s$
takes a physically allowed value, as is the case in the initial fixed–$t$
single dispersive integral of $\mathcal{M}$. Note that contributions from
heavy intermediate states are suppressed by: 1) the $1/s^{\prime}$ integral
kernel and 2) a reduced integration region.

In summary, from this theoretical framework it is possible to compute
scattering amplitudes starting from (31) and then obtaining the spectral
functions from experimental data for subprocesses via the optical theorem and
the Schwarz reflection principle999This sentence implicitly claims that the
analytic properties of an amplitude can be obtained entirely from its
dynamics, that is, from the intermediate states it allows. This claim is known
as the “Mandelstam hypothesis” and it is related to the causality requirement
of the theory..

####  1.3.3 Tensor decomposition of $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$

Dispersive integrals are suitable for scalar functions and
$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ is clearly not one. Hence, it is
necessary to decompose it, like we did for the electromagnetic vertex, into
form factors, which can then be dispersively represented. Nevertheless, the
Mandelstam hypothesis cannot be considered to apply in general for these
scalar coefficients. The key point behind this is that the tensor structures
of the decomposition may have kinematic singularities and/or zeroes, which
have to be cancelled by corresponding terms in the form factors, because the
amplitude does not have such terms according to the Mandelstam hypothesis. In
such cases, zeroes (singularities) change the asymptotic (analytic) behaviour
of the coefficients and this has an impact on the dispersion relation in the
form of subtraction constants in (31) and spurious residues. Such input has to
be determined experimentally and its presence in (31) further hinders the
computations101010It may even introduce ambiguities in the soft photon limit
of $\text{Im}\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ and its derivatives. See
[49].. In order to avoid these issues it is necessary to find a tensor
decomposition of the amplitude free of kinematic singularities and/or zeroes,
which is called a Bardeen–Tarrach–Tung (BTT) decomposition and was found
recently [31] for $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$. Let us review the main
steps followed in [31].

The only covariant objects which we can work with are the metric and the
momenta of the four photons. There are 138 possible combinations of said
objects. The most general structure is therefore:

$\begin{split}\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}&=g^{\mu_{1}\mu_{2}}g^{\mu_{3}\mu_{4}}\Pi^{1}+g^{\mu_{1}\mu_{3}}g^{\mu_{2}\mu_{4}}\Pi^{2}+g^{\mu_{1}\mu_{4}}g^{\mu_{3}\mu_{2}}\Pi^{3}\\\
&+\sum_{k,l}g^{\mu_{1}\mu_{2}}q_{k}^{\mu_{3}}q_{l}^{\mu_{4}}\Pi^{4}_{kl}+\sum_{j,l}g^{\mu_{1}\mu_{3}}q_{j}^{\mu_{2}}q_{l}^{\mu_{4}}\Pi^{5}_{jl}+\sum_{j,k}g^{\mu_{1}\mu_{4}}q_{k}^{\mu_{3}}q_{j}^{\mu_{2}}\Pi^{6}_{jk}\\\
&+\sum_{i,l}g^{\mu_{3}\mu_{2}}q_{i}^{\mu_{1}}q_{l}^{\mu_{4}}\Pi^{7}_{il}+\sum_{i,k}g^{\mu_{4}\mu_{2}}q_{k}^{\mu_{3}}q_{i}^{\mu_{1}}\Pi^{8}_{ik}+\sum_{i,j}g^{\mu_{3}\mu_{4}}q_{i}^{\mu_{1}}q_{j}^{\mu_{2}}\Pi^{9}_{ij}\\\
&+\sum_{i,j,k,l}q_{i}^{\mu_{1}}q_{j}^{\mu_{2}}q_{k}^{\mu_{3}}q_{l}^{\mu_{4}}\Pi^{10}_{ijkl}\;,\end{split}$
(36)

where the indices are $i\in\\{2,3,4\\}$, $j\in\\{1,3,4\\}$, $k\in\\{1,2,4\\}$
and $l\in\\{1,2,3\\}$111111Each index may be any three–element subset of the
four photons momenta. This choice in particular is due to [78] and is useful
to determine constraints from crossing symmetry.. There are no kinematic
singularities in the scalar coefficients untill this point. However there are
kinematic zeroes coming from two constraints that we have not yet explicitly
accounted for: gauge invariance and crossing symmetry.

Gauge invariance may be explicitly imposed by projecting each one of the
Lorentz indices of the amplitude onto the orthogonal space of the associated
virtual photon momentum. To this end the following projectors may be used
[79]:

$I_{12}^{\mu\nu}=g^{\mu\nu}-\frac{q_{1}^{\mu}q_{2}^{\nu}}{q_{1}\cdot
q_{2}},\hskip
28.45274ptI_{34}^{\mu\nu}=g^{\mu\nu}-\frac{q_{3}^{\mu}q_{4}^{\nu}}{q_{3}\cdot
q_{4}}\;,$ (37) $\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$
$\displaystyle=I_{12}^{\mu_{1}^{\prime}\mu_{1}}I_{12}^{\mu_{2}\mu^{\prime}_{2}}I_{34}^{\mu_{3}^{\prime}\mu_{3}}I_{34}^{\mu_{4}\mu^{\prime}_{4}}\Pi_{\mu^{\prime}_{1}\mu^{\prime}_{2}\mu^{\prime}_{3}\mu^{\prime}_{4}}\;,$
(38)

and only 43 tensor structures remain, which means that we have taken into
account 95 constraints coming from gauge invariance.

Of course these are not the only projectors suitable for the job. We could
have given each momentum its own projector. For example:

$\displaystyle
g^{\mu_{1}\mu_{2}}-\frac{q_{1}^{\mu_{1}}q_{1}^{\mu_{2}}}{q_{1}^{2}}\hskip
28.45274pt\text{or}\hskip
28.45274ptg^{\mu_{1}\mu_{2}}-\frac{q_{2}^{\mu_{1}}q_{2}^{\mu_{2}}}{q_{2}^{2}}\;.$

However, by applying any transverse projectors to
$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ we are introducing kinematic
singularities of the type $1/q_{i}\cdot q_{j}$. Thus, the less projectors we
introduce and the simpler they are, the less types of kinematic singularities
will be introduced.

In our case, two appearances of $I_{12}$ and two of $I_{34}$ in the tensor
structures lead to poles associated to all the possible combinations of
$q_{1}\cdot q_{2}$ and $q_{3}\cdot q_{4}$ with up to two repetitions of each.
This can be solved by building linear combinations of these singular tensor
structures such that the poles cancel. The precise procedure proposed by [79]
deals with the poles by decreasing singularity order. In the first step poles
of the form $(q_{1}\cdot q_{2})^{2}(q_{3}\cdot q_{4})^{2}$ are eliminated,
first by linear combinations and, once this is not possible, by multiplying
them by $q_{1}\cdot q_{2}$ or $q_{3}\cdot q_{4}$. The next step is to deal
with the single–double poles in the same way, that is, $(q_{1}\cdot
q_{2})(q_{3}\cdot q_{4})^{2}$ and $(q_{1}\cdot q_{2})^{2}(q_{3}\cdot q_{4})$.
Single–single poles then come and so forth until there are no kinematic
singularities left. The decomposition obtained after this procedure is
completed is not yet suitable for a dispersive representation because the
tensor “basis” found is actually not linearly independent in $q_{1}\cdot
q_{2}=0$ or $q_{3}\cdot q_{4}=0$ nor it actually spans the complete space of
possible gauge invariant tensors for the HLbL amplitude in those cases. This
is due to the existence of 11 linear combinations of the basis elements that
are proportional to $q_{1}\cdot q_{2}$ and/or $q_{3}\cdot q_{4}$ and some new
tensor structure, thus introducing the linear dependence and span issues just
mentioned.121212This phenomenon was first described by Tarrach [80] for the
BTT decomposition of the $\gamma\gamma\rightarrow\pi\pi$ process. These new 11
structures have to be added to the previous 43 to build a set that spans all
relevant gauge–invariant tensors even at $q_{1}\cdot q_{2}=0$ and $q_{3}\cdot
q_{4}=0$ [31]:

$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}=\sum_{i}\hat{T}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}_{i}\hat{\Pi}_{i}\,.$
(39)

Since the set is linearly dependent, the definition of the scalars
$\hat{\Pi}_{i}$ is obviously redundant.

It is not surprising that the 43–element basis does not work as expected in
the singular points of the projectors (37), but where do the new tensors come
from? To shed light on that issue let us consider, for example,
$\Pi^{7}_{21}$, which is one of the form factors that disappears after using
the gauge–invariance projectors. That means that it has been replaced by a sum
of other form factors using a constraint equation. It can be explicitly
obtained by studying the coefficient of $g^{\mu_{3}\mu_{2}}q_{1}^{\mu_{4}}$ in
the equation
$q_{1\mu^{\prime}_{1}}\Pi^{\mu^{\prime}_{1}\mu_{2}\mu_{3}\mu_{4}}=0$, which
yields the constraint:

$0=\Pi^{3}+(q_{1}\cdot q_{2})\Pi^{7}_{21}+(q_{1}\cdot
q_{3})\Pi^{7}_{31}+(q_{1}\cdot q_{4})\Pi^{7}_{41}\;.$ (40)

It is clear that $\Pi^{7}_{21}$ cannot be simplified using this constraint
when $q_{1}\cdot q_{2}=0$, thus the simplified tensor structures obtained
using this constraint no longer span the amplitude fully. Furthermore, the
remaining three form factors in (40) cease to be linearly independent when
$q_{1}\cdot q_{2}=0$, thus hinting at the other issue we found. One might
think that this is just an unfortunate result due to a poor solution of the
constraint system, i.e. bad choice of projectors, and that we should look for
one that solves for $\Pi^{3}$ instead of $\Pi_{21}^{7}$, for example.
Unfortunately, this is not possible: it is not hard to see that for all four
gauge invariance constraint equations only $\Pi^{1}$, $\Pi^{2}$ and $\Pi^{3}$
can be replaced without compromising the kinematics of the problem and there
are 95 different constraint equations, so there will always be at least 92
replacements that will have problematic limits.

Finally, it is worth noting that, thanks to the choice of the $i,j,k,l$ in
(36), the 54 tensor “basis”131313It is a basis in the sense that it spans the
tensor structures of the amplitude, but it is not linearly independent. is
closed under crossing. In fact, only 7 of the 54 tensors are actually
independent in terms of crossing transformations141414There are some elements
that are actually crossing antisymmetric, but the corresponding kinematic zero
does not affect the computation of $a_{\mu}$ [31]..

###  1.4 Master formula for the HLbL contribution to $a_{\mu}$

In this section we obtain $a_{\mu}^{\text{HLbL}}$ from the scalar coefficients
of (39), which is key to connect it to experimental data in the low energy
regime via dispersion relations.

Since the tensors $\hat{T}_{i}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ from (39) have
all the kinematic constrains and symmetries incorporated, they must vanish in
the soft photon limit like $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ does, as
argued previously. Besides, it is possible to find
$\hat{T}_{i}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ such that in the limit
$q_{4}\rightarrow 0$ the derivative of 35 of this tensors vanish, which leads
us from (28) to:

$\displaystyle a_{\mu}^{\text{HLbL}}$
$\displaystyle=e^{6}\int\frac{d^{4}q_{1}}{(2\pi)^{4}}\int\frac{d^{4}q_{2}}{(2\pi)^{4}}\frac{1}{q_{1}^{2}}\frac{1}{q_{2}^{2}}\frac{1}{(q_{1}+q_{2})^{2}}\frac{1}{(p+q_{1})^{2}-m^{2}}\frac{1}{(p-q_{2})^{2}-m^{2}}\,\sum_{i}^{19}\hat{T}_{i}\hat{\Pi}_{i}\;,$
$\displaystyle\hat{T}_{i}$
$\displaystyle\equiv\frac{1}{48m}\operatorname{Tr}\Big{\\{}(\not{p}+m)[\gamma_{\nu_{4}},\gamma_{\mu_{4}}](\not{p}+m)\gamma_{\mu_{1}}(\not{p}+\not{q}_{1}+m)\gamma_{\mu_{2}}(\not{p}-\not{q}_{2}+m)\gamma_{\mu_{2}}\Big{\\}}\frac{\partial}{\partial
q_{4\nu_{4}}}\hat{T}_{i}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}\Big{|}_{q_{4}\rightarrow
0}\;.$

The objects $\hat{T}_{i}$ act as kernels for the two loop integral. Their
number can be further reduced to 12 by harnessing the symmetry of the integral
and some of the kernels under the $q_{1}\leftrightarrow-q_{2}$ exchange, which
implies that some pairs of kernels actually give the same result and can be
absorbed into one.

At this point $a_{\mu}^{\text{HLbL}}$ seems to depend on $p$, but we know from
momentum conservation that $F_{2}=F_{2}(q_{4}^{2})$ and hence
$a_{\mu}^{\text{HLbL}}$ is just a number, not a function of any momentum. We
can remove the spurious dependence on $p$ through angular averages. Let us
start by performing a Wick rotation, in which we essentially render all
four–vectors’ time components imaginary. This causes all scalar products to
acquire a minus sign and renders the corresponding metric euclidean,
transforming the two loops accordingly. This has non–trivial consequences,
because it amounts to a rotation of the real axis of time–component
integration into the imaginary one. In case that the region swept by such
rotation contains singularities, the resulting contributions have to be taken
into account accordingly. However, for $a_{\mu}$ there no such issues [31] and
the Wick rotation may be performed without problems. The Wick–rotated version
of the momenta is represented by $Q_{1}$, $Q_{2}$ and $P$151515This capital
“P” still represents the initial momentum of the muon. It should not be
confused with $p+p^{\prime}$, as was represented in previous sections.,
$\hat{P}$ and $\hat{Q}_{i}$ represent their unit vectors, while $|P|$ and
$|Q_{i}|$ represent their norm. Then it is possible to remove the spurious
dependence on $P$ by averaging $a_{\mu}^{\text{HLbL}}$ over all possible
orientations of $P$:

$a_{\mu}^{\text{HLbL}}=\int\frac{d^{4}\Omega(P)}{2\pi^{2}}a_{\mu}^{\text{HLbL}}$
(41)

Wick rotated propagators strongly resemble the generating function of
Gegenbauer polynomials [81], which allows to represent each one of the former
as a linear combination of the latter. Finally, the integrals are found using
the polynomials’ orthogonality [82]. After performing the corresponding
average, it is possible to perform five of the six four–dimensional angular
integrals on $Q_{1}$ and $Q_{2}$ and the final result is:

$\displaystyle a_{\mu}^{\text{HLbL}}$
$\displaystyle=\frac{2\alpha^{3}}{3\pi^{2}}\int_{0}^{\infty}dQ_{1}\int_{0}^{\infty}dQ_{2}\int_{-1}^{1}d\tau\sqrt{1-\tau^{2}}|Q_{1}|^{3}|Q_{2}|^{3}\times\sum_{i}^{12}T_{i}\overline{\Pi}_{i}\;.$
(42)

There are three aspects of this last step that are worth noting:

  * •

The integral over $Q_{2}$ in spherical coordinates is considered in the first
place. It is possible to take any four–momentum as a reference for the angular
integral; it does not matter because the integrals will go over all the
possible values anyway. We take $Q_{1}$ as a reference.

  * •

The integrand is only dependent on one angle (in
$\tau=\hat{Q}_{1}\cdot\hat{Q}_{2}$) and it is therefore convenient to assign
$\tau$ as one of the three euclidean angles over which the angular integrals
of the four momentum $Q_{2}$ is performed. It is relevant which of these three
angle we are referring to because it will determine what the Jacobian will
look like. In the master formula there is a term $1-\tau^{2}$, a sine squared,
which means $\tau$ does not represent neither polar nor the azimuthal angle of
the three dimensional sphere embedded in the four dimensional space. Thus, the
angular integral on $Q_{2}$ yields:

$\displaystyle\int d\tau\sqrt{1-\tau^{2}}\int d\theta
d\phi\sin{\theta}=4\pi\int d\tau\sqrt{1-\tau^{2}}\;,$ (43)

where $\theta$ and $\phi$ represent the three–dimensional polar and azimuthal
angles of the four–dimensional $Q_{2}$ space.

  * •

Once the angular integrals on $Q_{2}$ have been performed there is no
dependence on $\tau$ or another angle left on the integrand. This means that
we can perform the four–dimensional solid angle integral on $Q_{1}$ which
yields $2\pi^{2}$.

##  2 Review of low energy contributions to $a_{\mu}^{\text{HLbL}}$

Up to this point we have presented an approach to compute low–energy
contributions to $a_{\mu}^{\text{HLbL}}$ based on a dispersive description of
the HLbL scattering amplitude. Now let us review the results obtained when
such approach is put to use.

As we stated in section 1.3, the dispersive calculation of an amplitude offers
a way to establish a hierarchy of contributions for the intermediate states
that enter the computation via the unitarity relation of (30). Heavier
intermediate states induce cuts that appear further to the right of the
dispersive integration region compared to their lighter counterparts and
therefore these former’s contributions are expected to be smaller. The
$1/(s-s^{\prime})$, $1/(t-t^{\prime})$ or $1/(u-u^{\prime})$ weights further
suppress such contributions. In the context of HLbL low energy scattering, the
most relevant intermediate states are expected to be the lightest one– or
two–hadron intermediate states.

###  2.1 One–particle intermediate states contribution to $a_{\mu}$

Baryons are not suitable single–particle intermediate states of the HLbL
scattering amplitude due to baryonic number conservation. Thus, they can only
contribute from two–particle intermediate states onwards. However, since the
proton, the lightest baryon, has a mass of $\sim 938.3$ MeV [83], its
corresponding two–particle state has a mass of about $1.9$ GeV, which is
already in the perturbative regime of QCD. This means that baryon
contributions can and should be accounted for perturbatively.

The lightest mesons with masses up to $\sim 1\leavevmode\nobreak\ \text{GeV}$
are $\pi^{0}$, $\pi^{\pm}$, $K^{\pm}$, $K^{0}$, $\eta$ and
$\eta^{\prime}$161616$\eta$ and $\eta^{\prime}$ are actually not stable in QCD
and therefore they do not appear in the unitarity relation (30) as a one-
particle intermediate state, however, its decay width is small enough to be
considered as a one.. Charge conservation obviously prohibits the appearance
of charged pions and kaons as one–particle intermediate states of HLbL.
Furthermore, a single neutral kaon cannot be produced since strong
interactions conserve strangeness. The contributions from the $\pi^{0}$,
$\eta$ and $\eta^{\prime}$ intermediate states to $a_{\mu}^{\text{HLbL}}$ are
of the form:

$\displaystyle a_{\mu}^{P-\text{Pole}}$
$\displaystyle=\frac{2\alpha^{3}}{3\pi^{2}}\int_{0}^{\infty}dQ_{1}\int_{0}^{\infty}dQ_{2}\int_{-1}^{1}d\tau\sqrt{1-\tau^{2}}|Q_{1}|^{3}|Q_{2}|^{3}\Big{(}T_{1}\overline{\Pi}_{1}^{P}+T_{2}\overline{\Pi}_{2}^{P}\Big{)}\;,$
(44) $\displaystyle\overline{\Pi}_{1}^{P}$
$\displaystyle=-\frac{\mathcal{F}_{P\gamma^{*}\gamma^{*}}(-Q_{1}^{2},-Q_{2}^{2})\mathcal{F}_{P\gamma^{*}\gamma^{*}}(-Q_{3}^{2},0)}{Q_{3}^{2}+M_{P}^{2}}\;,$
$\displaystyle\overline{\Pi}_{2}^{P}$
$\displaystyle=-\frac{\mathcal{F}_{P\gamma^{*}\gamma^{*}}(-Q_{1}^{2},-Q_{3}^{2})\mathcal{F}_{P\gamma^{*}\gamma^{*}}(-Q_{2}^{2},0)}{Q_{2}^{2}+M_{P}^{2}}\;,$

where $P$ states for the corresponding meson in the one–particle intermediate
state, $M_{P}$ is the meson’s mass and $\mathcal{F}_{P\gamma^{*}\gamma^{*}}$
represents its doubly virtual transition form factor, which is defined by:

$\displaystyle i\int d^{4}x\;e^{-iq_{1}x}\langle
0|T\\{J_{\mu}(x)J_{\nu}(0)\\}|P(q_{1}+q_{2})\rangle$
$\displaystyle=\epsilon_{\mu\nu\lambda\rho}q_{1}^{\lambda}q_{2}^{\rho}\mathcal{F}_{P\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})\;.$
(45)

The main hurdle for the computation of $a_{\mu}^{P-\text{pole}}$ is to obtain
reliable values for the transition form factor, specially for the double
virtual sector, for which there is usually little experimental data available
[84]. Furthermore, data available below $1$ GeV is particularly scarce, which
means that a fit extrapolation is required in order to obtain the transition
form factors at that range. This is specially troublesome because the main
contribution to (44) actually comes from the low energy region. For the
neutral pion, however, the negative effects of these issues have been
addressed by a dispersive reconstruction of the transition form factor [33,
85, 86] starting from (45). When $\mathcal{F}_{\pi^{0}\gamma^{*}\gamma^{*}}$
is computed in that fashion, the corresponding contribution to
$a_{\mu}^{\text{HLbL}}$ is:

$\displaystyle a_{\mu}^{\pi^{0}-\text{pole}}(\text{dispersive})$
$\displaystyle=63.0^{+2.7}_{-2.1}\times 10^{-11}\;.$ (46)

There are of course alternate approaches for evaluating
$\mathcal{F}_{\pi^{0}\gamma^{*}\gamma^{*}}$ such as the Canterbury
Approximants (CA) [87, 84, 88], which reproduces the low energy behaviour of
the transition form factor approximating it via rational functions of
polynomials in $Q_{1}^{2}$ and $Q_{2}^{2}\,$171717For the singly–virtual
transition form factor the Padé Approximants approach is used, which are
essentially the univariate version of the CA.. The CA computation of the pion
transition form factor yields the result [7]
$a_{\mu}^{\pi^{0}-\text{pole}}(\text{CA})=63.2(2.7)\times 10^{-11}$, which is
very much compatible with the dispersive one.

For $\eta$ and $\eta^{\prime}$ a dispersive computation of the corresponding
form factor in its doubly virtual region is not yet available, although the
framework has been established for the singly–virtual one [89]. Furthermore,
progress has been made in the doubly–virtual case as well [90, 91]. In any
case, until a full dispersive computation of the transition form factor is
obtained the CA approach offers an option whose reliability is supported by
the excellent compatibility with the dispersive results in the case of
$\pi^{0}$. For the $\eta$ and $\eta^{\prime}$ pole contributions the CA
approach concludes [84]:

$a_{\mu}^{\eta-\text{pole}}(\text{CA})=16.3(1.4)\times 10^{-11}\;,\hskip
28.45274pta_{\mu}^{\eta^{\prime}-\text{pole}}(\text{CA})=14.5(1.9)\times
10^{-11}\;.$ (47)

It is clear that the neutral pion contribution is by far the largest from the
one–particle intermediate states. Since the $\pi^{0}$ is lighter than the
$\eta$ and it in turn is lighter than the $\eta^{\prime}$, these results for
the contribution of each of these intermediate states agree with the expected
hierarchy mentioned at the beginning of this section. In summary, the
one–particle intermediate states contribution to $a_{\mu}$ is [7]:

$\displaystyle a_{\mu}^{\pi^{0}+\eta+\eta^{\prime}}$
$\displaystyle=93.8^{+4.0}_{-3.6}\times 10^{-11}\;.$ (48)

###  2.2 Two–particle intermediate states contribution to $a_{\mu}$

The lightest two–meson intermediate states with mass up to $\sim 1$ GeV from
lightest to heaviest are $\pi^{0}\pi^{0}$, $\pi^{+}\pi^{-}$, $\eta\pi^{0}$,
$K^{+}K^{-}$ and $K^{0}\overline{K}^{0}$. Note that, again due to strangeness
conservation in QCD, no intermediate states with a single kaon is allowed.

The dispersive frame work for the computation of the two–pions contribution
has been studied in [31]. In there, by performing the unitarity cut across the
s–channel, this amplitude is split into two $\gamma\gamma\rightarrow\pi\pi$;
one with two virtual photons and one with one virtual and one real soft
photon. Then, this two amplitudes are split again by performing the unitarity
cut across the t–channel, which is the application of a second dispersion
relation that leads to the Mandelstam representation that was described in
section 1.3. Although in this second cut there is the possibility for
multi–pion intermediate states to appear, the biggest contribution is expected
to come from the pion pole (lightest intermediate state) as well. This
introduces the pion–box topology as an intermediate state.181818It is worth
noting that due to angular momentum conservation and Bose symmetry, a photon
does not couple to two neutral pions and therefore the
$\gamma^{*}\gamma^{*}\rightarrow\pi^{0}\pi^{0}$ does not allow a neutral pion
intermediate state across the t–channel. Since the four pieces in which the
original HLbL is split are essentially pion electromagnetic vertices, it is
reasonable to expect the pion–box topology contribution to $a_{\mu}$ to be
proportional to four pion electromagnetic vector form factors
$\mathcal{F}^{V}_{\pi}(q_{i}^{2})$. Furthermore in [31] this relation is taken
further and via explicit computation of the double spectral functions
associated to the pion–box contribution it was shown that the pion–box
contribution to $\Pi^{\mu\nu\lambda\rho}$ is in fact equal to the scalar QED
pion loop amplitude multiplied by four pion vector form factors due to each of
the four off-shell photons. The mathematical expression of these statements
is:

$\displaystyle\Pi_{i}^{\pi-\text{box}}$
$\displaystyle=\mathcal{F}^{V}_{\pi}(q_{1}^{2})\mathcal{F}^{V}_{\pi}(q_{2}^{2})\mathcal{F}^{V}_{\pi}(q_{3}^{2})\mathcal{F}^{V}_{\pi}(q_{4}^{2})\times\frac{1}{\pi^{2}}\Bigg{(}\int_{c_{1}}^{\infty}ds^{\prime}\int_{c_{2}(s^{\prime})}^{\infty}dt^{\prime}\frac{\Delta_{t^{\prime}}\Delta_{s^{\prime}}\mathcal{M}(s^{\prime},t^{\prime})}{(s^{\prime}-s)(t^{\prime}-t)}$
$\displaystyle+\int_{c_{1}}^{\infty}ds^{\prime}\int_{c_{2}(s^{\prime})}^{\infty}du^{\prime}\frac{\Delta_{u^{\prime}}\Delta_{s^{\prime}}\mathcal{M}(s^{\prime},u^{\prime})}{(s^{\prime}-s)(u^{\prime}-u)}+\int_{c_{1}}^{\infty}du^{\prime}\int_{c_{2}(u^{\prime})}^{\infty}dt^{\prime}\frac{\Delta_{t^{\prime}}\Delta_{s^{\prime}}\mathcal{M}(t^{\prime},u^{\prime})}{(u^{\prime}-u)(t^{\prime}-t)}\Bigg{)}\;,$
(49)

where $\mathcal{M}$ is the scalar QED one loop amplitude with pions. The
corresponding double spectral functions can of course be obtained used
Cutkosky’s rules to perform unitarity cuts to loops, however, it is much more
efficient to compute the light by light scattering amplitude directly in
perturbative scalar QED and then insert into the master formula (42) with the
appropriate vector form factors. This allows for the contribution of the pion
box to be written in terms of compact and well known Feynman parameters
integrals. In conclusion, the pion box contribution to $a_{\mu}^{\text{HLbL}}$
reads:

$\displaystyle a_{\mu}^{\pi-\text{box}}$
$\displaystyle=\frac{2\alpha^{3}}{3\pi^{2}}\int_{0}^{\infty}dQ_{1}\int_{0}^{\infty}dQ_{2}\int_{-1}^{1}d\tau\sqrt{1-\tau^{2}}|Q_{1}|^{3}|Q_{2}|^{3}\mathcal{F}^{V}_{\pi}(-Q_{1}^{2})\mathcal{F}^{V}_{\pi}(-Q_{2}^{2})\mathcal{F}^{V}_{\pi}(-Q_{3}^{2})\sum_{i}^{12}T_{i}\overline{\Pi}_{i}^{sQED}$
$\displaystyle=-15.9(2)\times 10^{-11}\;,$ (50)

where $\overline{\Pi}_{i}^{sQED}$’s are the scalar coefficients of the BTT
decomposition of the scalar QED one–loop light–by–light scattering amplitude,
which obviously is a subset of the HLbL BTT structures. In the previous
equation only three vector form factors appear, instead of four, because one
of the photons of this HLbL process is actually on–shell.

Regarding intermediate states heavier than one pion in the t–channel cut of
the $\gamma\gamma\rightarrow\pi\pi$ subprocess, they can be classified into
two categories: one with one pion–pole contribution in one subamplitude and a
multiparticle cut in the other, and one where both subprocesses contain
multiparticle cuts. The computation of this contributions can be performed
through a partial wave expansion of the amplitude. One of the advantages of
this approach is that the unitarity relation (30) is diagonal in the helicity
partial wave basis, that is, helicity partial waves of HLbL
$h^{J}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}$ are only connected by
unitarity with the $\gamma^{*}\gamma^{*}\rightarrow\pi\pi$ partial waves
$h^{J}_{\lambda_{i}\lambda_{j}}$ with the same total angular momentum $J$ and
photon helicities $\lambda_{i}$:

$\begin{split}\text{Im}^{\pi\pi}\;h^{J}_{\lambda_{1},\lambda_{2},\lambda_{3},\lambda_{4}}=\frac{\sigma_{\pi}}{16\pi
S}h^{J}_{\lambda_{1}\lambda_{2}}h^{*J}_{\lambda_{3}\lambda_{4}}\;,&\hskip
14.22636ptH_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}=\sum_{J}(2J+1)d_{m_{1}m_{2}}^{J}(z)h^{J}_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}\;,\\\
H_{\lambda_{i}\lambda_{j}}=\sum_{J}(2J+&1)d_{m0}^{J}(z)h^{J}_{\lambda_{i}\lambda_{j}}\;,\end{split}$
(51) $\begin{split}m=|\lambda_{i}-\lambda_{j}|\;,&\hskip
28.45274ptm_{1}=|\lambda_{1}-\lambda_{2}|\;,\\\
m_{2}=|\lambda_{3}-\lambda_{4}|\;,&\hskip
28.45274pt\sigma_{\pi}=\sqrt{1-\frac{4M^{2}_{\pi}}{s}}\;,\end{split}$

where $H_{\lambda_{1}\lambda_{2}\lambda_{3}\lambda_{4}}$ and
$H_{\lambda_{i}\lambda_{j}}$ are helicity amplitudes for the HLbL and the
$\gamma^{*}\gamma^{*}\rightarrow\pi\pi$ processes, respectively, which are
obtained by contracting polarization vectors of a given helicity with the
corresponding process amplitude. Furthermore, $d_{m_{1}m_{2}}^{J}$ are
Wigner’s matrices, $z$ is the scattering angle and $M_{\pi}$ represents the
pion’s mass. Helicity partial waves are, however, not suitable for a double
Mandelstam representation, because they contain kinematic singularities (see
section 3 of chapter 7 in [92]). To retain the advantages of the partial wave
expansion, but solve the kinematic singularities issues it is necessary to
relate them to the BTT basis described in section 1.3. This is done by first
inverting the linear relation that is obtained between helicity amplitudes and
BTT scalar functions via its definition from contraction with photon
polarization vectors. Then helicity partial waves are projected out of the
corresponding helicity amplitude by using orthogonality properties of Wigner’s
matrices. As we mentioned in section 1.3 [77], the BTT decomposition of
$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ actually trades linear independence in
order to span the whole amplitude, therefore introducing redundancies. This
redundancies are constrained in the soft limit ($q_{4}\rightarrow 0$) by using
dispersive sum rules derived from the asymptotic properties of the scalar
coefficients [32]. Finally, the S–wave contribution to the two pion HLbL
intermediate states with one pion as and intermediate state in one
$\gamma^{*}\gamma^{*}\rightarrow\pi\pi$ subprocess and a multiparticle cut in
the other is [32]:

$\displaystyle a_{\mu}^{\pi\pi,\;\pi-\text{pole LHC}}=-8(1)\times 10^{-11}\;,$
(52)

where “LHC” stands for “left hand cut”. Contributions from higher angular
momentum partial waves is difficult due to the lack of experimental input for
the doubly virtual subprocess $\gamma^{*}\gamma^{*}\rightarrow\pi\pi$.
Furthermore, heavier singularities in the left hand cut have to be taken into
account for the computation of higher angular momentum partial waves, but in
such case some assumptions that simplified the computation of (52) are no
longer valid. Therefore, the dispersive framework developed in [32] would have
to be studied from a more general perspective.

In [7] it is proposed that the comparison between the size of the
$\gamma\gamma\rightarrow\pi\pi$ cross section and $\gamma\gamma\rightarrow MM$
cross section for some given processes gives input to determine which
two–particle intermediate state contribution $MM$ may be relevant beyond the
$\pi\pi$ contribution. Following this approach, the relevant heavier
intermediate states are $\pi^{0}\eta$ and $K^{+}K^{-}$, which are actually the
next two–particle intermediate states in terms of mass. For the two kaon
contribution, the same arguments of the pion box apply and therefore, via the
scalar QED kaon loop and using vector meson dominance to obtain the
electromagnetic form factor, the corresponding contribution has bee computed
[7]:

$\displaystyle a_{\mu}^{KK\text{--loop}}$ $\displaystyle=-0.5(1)\times
10^{-11}\;.$ (53)

For the $\pi^{0}\eta$ the real version of the process a dispersive framework
has been developed [93] and the results performed well against data of the
$\eta\rightarrow\pi^{0}\gamma\gamma$ crossed process. Work on the singly and
doubly virtual processes is still under way [94].

###  2.3 Intermediate states with more than two particles and heavier than
$KK$

The dispersive computation of contributions heavier than the ones we presented
previously, that is, with masses between $1$ GeV and $2$ GeV, is more
difficult due to the more complex unitarity diagrams than come into play.
Perhaps the most straightforward upgrade is the contribution from D–191919Each
pion has isospin equal to one, therefore the possible values for total isospin
are $|1-1|=0$ and $1+1=2$. This means that only even partial waves are
allowed. and higher waves from the two pion contribution, for which there is a
pretty much developed framework [34]. Contributions from states heavier that
two pions are expected to be small. For example, the $K^{+}K^{-}$ intermediate
state has mass $\sim 1$ GeV and its contribution to $a_{\mu}$ is already below
$10^{-11}$. Consequently, heavier intermediate states contributions are
expected to be non–negligible only if they are associated to resonant helicity
partial waves. To estimate such contributions a narrow–width resonance
approximation is followed in [95] (updated in [96]) and [97], which
essentially means that narrow resonances are treated as dispersive poles. In
summary, the contribution considering the $f_{0}(980\text{ MeV})$,
$a_{0}(980\text{ MeV})$, $f_{0}(1370\text{ MeV})$, $a_{0}(1450\text{ MeV})$,
$f_{0}(1500\text{ MeV})$, $f_{2}(1270\text{ MeV})$, $f_{2}(1565\text{ MeV})$,
$a_{2}(1320\text{ MeV})$ and $a_{2}(1700\text{ MeV})$ 202020In [95] and its
updated version [96] the first three (scalar) resonances and the last four
(tensor) resonances are considered. In [97] considers all the scalar
resonances mentioned. is:

$\displaystyle a_{\mu}^{\text{scalars}+\text{tensor}}$
$\displaystyle=-1(3)\times 10^{-11}\;,$ (54)

which introduces a significant uncertainty due to the unknown error introduced
by the narrow–width resonance approximation. In principle, such approximation
may be tested in the lightest scalar and tensor resonances against the
$\pi\pi$ and $KK$ contributions described previously. However, only a full
tower of resonant poles satisfy the sum rules found in [32], which is a
necessary condition for a tensor–basis independent contribution.

Axial vector mesons (meson bound states with total angular momentum equal to
1, but even parity) have not been considered in the previous paragraphs,
because212121Note that the Landau-Yang theorem only forbids the coupling of
two real photons to axial vector mesons, therefore in virtual HLbL it is
allowed. the lightest associated resonances are actually heavy compared with
other states: $a_{1}(1260\text{ MeV})$, $f_{1}(1285\text{ MeV})$ and
$f_{1}(1420\text{ MeV})$. They contribute to $a_{\mu}$ as resonant poles and
their respective transition form factors are obtained from experimental data
using a non–relativistic quark model computation [98, 99, 100, 101] which
simplifies the tensor structure of the $\gamma^{*}\gamma^{*}\rightarrow A$
amplitude, as described in [95]. The current agreed contribution (see [7]) is
obtained from the mean of three results [95, 97, 102] that consider the three
lightest axial vector mesons:

$\displaystyle a_{\mu}^{\text{axial vectors}}$ $\displaystyle=6(6)\times
10^{-11}\;.$ (55)

In summary, we have presented the basics of the computation of $a_{\mu}$ and
in particular $a_{\mu}^{\text{HLbL}}$, which led to the (42) master formula,
which is based on a BTT tensor decomposition of the HLbL amplitude with 54
elements such that its scalar coefficients are free of kinematic singularities
and zeroes. The tensor “basis” is actually not linearly independent in the
whole phase space due to 11 structures that need to be added by hand in order
to span the amplitude at some kinematic points $q_{1}\cdot q_{2}=0$ or
$q_{3}\cdot q_{4}=0$. Furthermore, in four space–time dimensions there are two
additional linear relations among the 138 tensors of (36) [103] which the 54
BTT structures inherite. Of course, the linear dependence of the set of
structures introduces redundancies in the definition of the associated scalar
coefficients. However, in [34] it was shown that these redundancies do not
affect the observables, as expected, due to a series of dispersive sum rules
obeyed by the form factors.

We also discussed contributions of intermediate states with masses up to $\sim
2$ GeV, which are expected to give the man contribution to
$a_{\mu}^{\text{HLbL}}$. However, it is as well necessary to understand the
behaviour of the HLbL tensor beyond that threshold. The high energy regime of
the HLbL contribution has been already reached from the light intermediate
states when the high energy part of the integral in (42) has been carried out
for high virtuality values of the transition and vector form factors. Although
these penetrations into the high energy regime are expected to take into
account most of the contribution, in this regime heavier states with more
complex topologies that are not taken into account dispersively can also
contribute. It is those contributions and other uses for computations in the
high energy regime of the QCD in the context of $a_{\mu}$ with which the next
sections will deal.

##  3 Operator product expansion of $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$

Computations in the high energy regime of QCD play several important roles in
the determination of $a_{\mu}$. First, for data–driven computations they
provide the asymptotic behaviour of transition and electromagnetic form
factors for the high energy “tail” of the integral in the master formula (42)
and dispersive integrals. Examples of these short distance constraints (SDC)
can be found in [104, 105, 106, 107, 108, 109] and in [110, 111] for the
specific case of the pion. Second, SDC can be found for the HLbL amplitude
itself, which are used to evaluate how much of this asymptotic behaviour is
recovered by dispersive computations and thus assess the uncertainty produced
by the missing high energy contributions. These come from the heavy
intermediate states with topologically more complex unitarity diagrams that
are ignored in the dispersive approach.

In this work we are interested in the SDC for the HLbL tensor. There are two
high energy regimes for the HLbL process with one real soft photon: one where
all three Euclidean virtualities are similarly high ($Q^{2}_{1}\sim
Q^{2}_{2}\sim Q^{2}_{3}\gg\Lambda^{2}_{QCD}$), (we refer to the Euclidean
virtualities $q_{i}^{2}=-Q_{i}^{2}$) and one where two are similarly high and
much greater than the third Euclidean virtuality ($Q^{2}_{1}\sim Q^{2}_{2}\gg
Q^{2}_{3}\sim\Lambda^{2}_{QCD}$ and crossed versions). It is worth to remind
that we refer to the Euclidean virtualities, because the considered photons
are far off-shell since when the space-time separation goes to zero the
corresponding momenta $q_{i}$ become space-like. Each of these regimes impose
asymptotic behaviour constraints on different subsets of BTT scalar
coefficients and therefore they allow for the independent evaluation of the
different sets of intermediate states. In this section we will focus on the
$Q^{2}_{1}\sim Q^{2}_{2}\sim Q^{2}_{3}\gg\Lambda^{2}_{QCD}$ regime of the HLbL
tensor, which we will study by performing an Operator Product Expansion (OPE)
where the soft photon is introduced as a background electromagnetic field, as
was done in [56].

###  3.1 OPE of $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$

A product of operators carrying very high momenta can be expressed in terms of
a linear combination of local operators carrying zero momentum with singular
coefficients that carry the momentum dependence of the original operator. The
mathematical expression of this statement for the product of two operators
$O_{1}$ and $O_{2}$ is:

$\displaystyle\int d^{4}xe^{-iqx}O_{1}(x)O_{2}(0)$
$\displaystyle=\sum_{n}C_{n}^{12}(q^{2})O_{n}(0)\;,$ (56)

where $C_{n}^{12}$ is a c–number function that decreases for large $Q^{2}$ and
is known as “Wilson coefficient”. The operator basis $\\{O_{n}\\}$ of the
expansion only admits elements with the same quantum numbers and symmetries of
the original operator product. The generalization for higher numbers of
operators with arbitrary tensor structures is straightforward. Such identity
is called operator product expansion and was originally introduced by Kenneth
Wilson [112]. By engineering dimensional analysis it can be noticed that the
Wilson coefficients vanish increasingly quicker for $O_{n}$ of higher mass
dimensions, which implies that the OPE has an implicit and useful hierarchy of
contributions in which simpler operators (with less derivatives and fields)
give bigger contributions than more complex operators. Since the operators on
both sides of (56) are renormalized at some scale $\mu$, a more careful study
finds that dimensional analysis must be done considering the anomalous
dimensions of the operators $O_{n}$, but for asymptotically free theories such
as QCD, at high energy the simple power counting suffices. It is also worth
noting that the relation in (56) is given in terms of operators, that is, it
does not depend on specific matrix elements an therefore Wilson coefficients
do not either. For further details see chapter 20 of volume II in [113] and
chapter 18 in [75].

To perform the OPE in (56) it is first necessary to fix the maximum number of
mass dimensions of the operators that are going to be considered and then the
task is to find all compatible operators with the same quantum numbers and
symmetries of the original operator product. Once the complete set of relevant
operators $O_{n}$ is known, each Wilson coefficient is found by choosing an
appropriate matrix element of the original operator product $O_{1}O_{2}$ and
then transforming it into the corresponding element of the operator basis
$O_{n}$. Such transformation is done by expanding the correlation function and
leaving the elements that form $O_{n}$ uncontracted. The result of the
contraction of the rest of the parts of the original product and the fields in
the correlation function vertices constitutes the Wilson coefficient of the
element $O_{n}$ for the operator basis.

From the definition of the OPE it is evident that it constitutes a very well
suited framework for the evaluation of the HLbL tensor in the $Q_{1}^{2}\sim
Q_{2}^{2}\sim Q_{3}^{2}\gg\Lambda_{QCD}^{2}$ regime. However, the limit
$q_{4}\rightarrow 0$ does not allow to include the fourth current of (24) in
the OPE. One could in principle start the construction of such OPE for the
four currents of $\Pi_{\text{HLbL}}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$, but it
does not work for our particular problem. For example, the Wilson coefficient
associated to the identity operator is just the HLbL tensor in perturbative
QCD, which involves an expansion in terms of the strong coupling constant and
the usual large logarithms for increasing powers:
$\alpha^{n}_{s}(\mu)\ln^{n}{\big{\\{}Q_{4}^{2}/\mu^{2}\big{\\}}}$, where $\mu$
represents the renormalization subtraction point in the $\overline{MS}$
scheme. In order for the logarithms not to blow up it is necessary to have
$\mu\sim Q_{4}$, but in such case $\alpha_{s}$ would enter the
non–perturbative domain of QCD and the expansion would be spoiled anyway.
Wilson coefficients for higher dimensional operators also suffer from infrared
singularities since they depend upon the $1/q_{4}^{2}$ propagator. Including
the fourth current $J^{\mu_{4}}$ into the OPE even though its momentum is not
large is not the only problem. Even if the OPE were performed only for the
three currents with high momenta, there would still be matrix elements of the
type $\langle 0|O_{n}J^{\mu_{4}}(q_{4})|0\rangle$, which cannot be
perturbatively computed in QCD. In general, these issues are different
consequences of the fact that perturbative QCD is not the correct framework to
describe the soft interaction that is required by $a_{\mu}$. It is therefore
necessary to perform the OPE of the three high–momentum currents only and also
take the fourth one into account properly. This can be done by letting the
soft photon be introduced by an external electromagnetic field instead of a
quark current and we will present such procedure in the next section of this
review222222This approach was first used in [114] in the context of the
computation of the magnetic moment of nucleons, then it was used in [115] for
the hadronic corrections to the electroweak contribution to $a_{\mu}$ and
finally it was again picked up in [55, 56, 57] for the HLbL tensor. A
pedagogical review of the framework is presented in [116]..

###  3.2 OPE of $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ in an electromagnetic
background field: A first look

A new object suitable for the OPE described at the end of the previous part
is:

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}}$ $\displaystyle=\frac{1}{e}\int
d^{4}x\int d^{4}y\;e^{-i(q_{1}x+q_{2}y)}\langle
0|T\,J^{\mu_{1}}(x)J^{\mu_{2}}(y)J^{\mu_{3}}(0)|\gamma(q_{4})\rangle=-\epsilon_{\mu_{4}}(q_{4})\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}\;,$

where the soft photon $q_{4}\rightarrow 0$ is included implicitly in the
initial state. In addition, this time $J^{\mu}$ makes reference to the
electromagnetic current of the three lightest quarks, namely, up, down and
strange or $u$, $d$ and $s$ and thus:

$\begin{split}J^{\mu}&=\overline{\Psi}\hat{Q}\gamma^{\mu}\Psi\hskip
56.9055pt\hat{Q}=\text{diag}\Big{(}\frac{2}{3},-\frac{1}{3},-\frac{1}{3}\Big{)}\;,\end{split}$
(57)

where $\hat{Q}$ is the charge matrix and now $\Psi$ is a vector of bispinors
with quark flavor and color indices, which are summed upon and suppressed in
the current.

Since $\langle 0|...|\gamma(q_{4})\rangle$ is an on–shell matrix element only
gauge invariant operators contribute to the OPE. From these only $F_{\mu\nu}$,
the field–strength tensor, contributes to first order in $q_{4}$, hence only
operators that have the same quantum numbers and symmetries of $F_{\mu\nu}$
are relevant. In summary, this means that at the first order in the external
electromagnetic field we have for the regime of high virtualities at hand:

$\begin{split}\Pi^{\mu_{1}\mu_{2}\mu_{3}}&\equiv
i\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\mu_{5}}_{F}(q_{1},q_{2})\langle
0|F_{\mu_{4}\mu_{5}}(0)|\gamma(q_{4})\rangle=q_{4\mu_{4}}\epsilon_{\mu_{5}}(q_{4})\Pi^{\mu_{1}\mu_{2}\mu_{3}[\mu_{4}\mu_{5}]}_{F}\;.\\\
\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}}{\partial
q_{4\mu_{5}}}\Bigg{|}_{q_{4}\rightarrow
0}&=\Pi^{\mu_{1}\mu_{2}\mu_{3}[\mu_{4}\mu_{5}]}_{F}\;.\end{split}$ (58)

We see here a confirmation of the antisymmetric nature of
$\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ arising from gauge invariance. We
can also see that the real object of interest for $a_{\mu}^{\text{HLbL}}$ is
actually $\Pi^{\mu_{1}\mu_{2}\mu_{3}[\mu_{4}\mu_{5}]}_{F}$, which is made from
the Wilson coefficients of the OPE. Note that t is explicitly free of any
$q_{4}$ dependence and therefore does no suffer from singularities at
$q_{4}^{2}\rightarrow 0$.

We already discussed that the OPE elements for $\Pi^{\mu_{1}\mu_{2}\mu_{3}}$
have the same Lorentz structure and symmetries of $F_{\mu_{1}\mu_{2}}$, that
is:

  * •

second rank antisymmetric tensor;

  * •

odd charge–conjugation parity (remember in this regard the famous Furry’s
theorem).

In [56] operators with these features and mass dimension up to $6$ are taken
into account and the rest are neglected. This choice is ultimately supported
by the fact the contribution of higher dimensional operators turns out to be
at least two orders of magnitude smaller than the leading order. It is however
also true that the non perturbative matrix elements of the dimension seven
operators are less known. We focus here, for simplicity of reading at this
step of our analysis, to the case of only one flavour and therefore on the
following list of operators

$\begin{split}S_{1,\mu\nu}&\equiv ee_{f}F_{\mu\nu}\;,\hskip
46.23572ptS_{2,\mu\nu}\equiv\overline{\Psi}\sigma_{\mu\nu}\Psi\;,\\\
S_{3,\mu\nu}&\equiv ig_{S}\overline{\Psi}G_{\mu\nu}\Psi\;,\hskip
28.45274ptS_{4,\mu\nu}\equiv
ig_{S}\overline{\Psi}\,\overline{G}_{\mu\nu}\gamma_{5}\Psi\;,\\\
S_{5,\mu\nu}&\equiv\overline{\Psi}\Psi\;ee_{f}F_{\mu\nu}\;,\hskip
28.45274ptS_{6,\mu\nu}\equiv\frac{\alpha_{s}}{\pi}G^{\alpha\beta}_{a}G^{a}_{\alpha\beta}\;ee_{f}F_{\mu\nu}\;,\\\
S_{7,\mu\nu}&\equiv
g_{S}\overline{\Psi}(G_{\mu\lambda}D_{\nu}+D_{\nu}G_{\mu\lambda})\gamma^{\lambda}\Psi+g_{S}\overline{\Psi}(G_{\nu\lambda}D_{\mu}+D_{\mu}G_{\nu\lambda})\gamma^{\lambda}\Psi\;,\\\
S_{\\{8\\},\mu\nu}&\equiv\alpha_{s}(\overline{\Psi}\Gamma\Psi\overline{\Psi}\Gamma\Psi)_{\mu\nu}\;,\end{split}$
(59)

where $\Psi$ represents again a quark field in a given flavour of electric
charge $ee_{f}$, the colour indices indices are implicitly summed upon,
$\Gamma$ represents a combination of Dirac gamma matrices, $D_{\nu}$
represents the gauge–covariant derivative, $G_{\mu\nu}^{a}$ represents the
gluon field strength tensor, $G_{\mu\nu}\equiv it^{a}G_{\mu\nu}^{a}$ and
$\overline{G}_{\mu\nu}\equiv\frac{i}{2}\epsilon^{\mu\nu\alpha\beta}G_{\alpha\beta}$.
Since the largest non–perturbative QCD energy scale is the perturbative
threshold $\Lambda_{QCD}$, then the contributions from operators with mass
dimension $d$ are expected to be suppressed like
$\Big{(}\frac{\Lambda_{QCD}}{Q_{i}}\Big{)}^{d-1}$232323Note that the mass
dimensions of the matrix element $\langle 0|...|\gamma\rangle$ of an operator
with mass dimension $d$ actually has mass dimension $d-1$., which means that
an $O\Big{(}\frac{\Lambda_{QCD}^{d}}{Q_{i}^{d}}\Big{)}$ error is introduced
when a cut–off dimension $d$ is imposed. When included in the high–energy
integration region of the master formula (42), the integration domain should
be used coherently in agreement with the mass dimension cut–off of the OPE.

Due to the quantum numbers of the $S_{i,\mu\nu}$ operators we will be able to
factorize the plane wave of the external electromagnetic field times the
non–zero expectation value in the true QCD vacuum of a Lorentz, gauge
invariant and charge-conjugation even operator, so called a condensate
$X_{i}^{S}$:

$\displaystyle\langle 0|S_{i,\mu\nu}|\gamma\rangle\equiv X_{i}^{S}\langle
0|F_{\mu\nu}|\gamma\rangle\;,$ (60)

which render the connection between the $S_{i\mu\nu}$ operators and $\Pi_{F}$
more evident.

Now that the operator basis is known up to dimension six, the next step is to
obtain the Wilson coefficients of its elements. Since these are local
operators, evaluated at $x=0$, the obvious step is to perform a Taylor
expansion of the field variables. However, such an expansion hinders the
computation of the Wilson coefficients because it involves terms that are not
even gauge covariant while the $S_{i,\mu\nu}$ are gauge invariant. Hence the
computation of the Wilson coefficients becomes very complex, particularly for
the gluon operators, as can be seen in [117, 118].

Instead, if one chooses the Fock–Schwinger gauge [119] for the photon and
gluon fields, defined by the constraint $A_{\mu}(x-x_{0})^{\mu}=0$ for some
constant $x_{0}$, the following expansions hold [120]:

$\begin{split}\Psi(x)&=\Psi(0)+x^{\mu_{1}}D_{\mu_{1}}\Psi(0)+\frac{1}{2!}x^{\mu_{1}}x^{\mu_{2}}D_{\mu_{1}}D_{\mu_{2}}\Psi(0)+...\;,\\\
A^{a}_{\alpha}(x)&=\frac{1}{2\times
0!}x^{\mu_{1}}G^{a}_{\mu_{1}\alpha}(0)+\frac{1}{3\times
1!}x^{\mu_{1}}x^{\mu_{2}}D_{\mu_{1}}G_{\mu_{2}\alpha}^{a}(0)+...\;,\\\
A_{\alpha}(x)&=\frac{1}{2\times
0!}x^{\mu_{1}}F_{\mu_{1}\alpha}(0)+\frac{1}{3\times
1!}x^{\mu_{1}}x^{\mu_{2}}D_{\mu_{1}}F_{\mu_{2}\alpha}(0)+...\;,\\\
\end{split}$ (61)

where $A^{a}$ represents the gluon field and we have chosen $x_{0}=0$ for
simplicity. The gauge–covariant form of these expansions offers a clear
advantage for the computation of Wilson coefficients associated to
gauge–invariant operators. On the other hand, the gauge–fixing constraint
evidently breaks traslational invariance and the propagator, being gauge
dependent, inherits this feature, which results in a rather complex expression
[121]. As a result the gain in simplicity due to explicit gauge covariance may
be lost in perturbative corrections due to the complexity of the propagators.

A way around this issue is to split the gauge field into a background part
that acts as a classical field that parameterizes the effects of the
non–perturbative QCD vacuum on perturbative dynamics dictated by the
asymptotic freedom and a fluctuation part that represents the perturbative
oscillations around the vacuum solution. The former leads to the non–zero
values $\langle S_{i,\mu\nu}\rangle$ and can be fixed by the radial gauge
constraint, but it is not a dynamical field variable and hence does not have
an associated propagator, and viceversa for the fluctuation part, which is
fixed by a different gauge constraint that we will explore in the next
section. Such a split can be represented as:

$A_{\mu}(x)=a_{\mu}(x)+A^{{}^{\prime}}_{\mu}(x)\;,\hskip
19.91684ptA^{a}_{\mu}(x)=a^{a}_{\mu}(x)+A^{{}^{\prime}a}_{\mu}(x)\;,$ (62)

where the unprimed variables are the classical fields an the primed ones
represent the fluctuations around this classical value. It is not needed to
include the specific form of the QCD vacuum fields but just to parameterize
them as external fields, as it is for instance done in the Coleman-Weinberg
approach. In the next section we present in detail the theoretical framework
that supports the background field method for these split fields.

###  3.3 OPE of $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ in an electromagnetic
background field: Theoretical framework

First, we will relate the results obtained with the split fields to the ones
obtained with the conventional approach. Let $Z[J_{n}]$ represent the
generating functional for the HLbL interaction:

$Z[J_{n}]=\int\prod_{n}\mathcal{D}\Phi_{n}\exp{i\Big{(}S[\Phi_{n}]+J_{n}\Phi_{n}\Big{)}}\;,$
(63)

where $\Phi$ represents different types of fields (fermion, vector, etc.) and
$n$ is a collective index that is meant to label all (color, flavor and
Lorentz) degrees of freedom. Complementary one can define a background
generating functional:

$\mathcal{Z}[J_{n},\phi_{n}]=\int\prod_{n}\mathcal{D}\phi^{\prime}_{n}\exp{i\Big{(}S[\phi_{n}+\phi^{\prime}_{n}]+J_{n}\phi^{\prime}_{n}\Big{)}}\;,$
(64)

where primed and unprimed variables represent quantum fluctuations and
background fixed fields, respectively, just as in the previous section. From
these objects one can obtain the corresponding quantum effective action by
performing a Legendre transform on the generator of connected diagrams:

$\displaystyle W[J_{n}]$ $\displaystyle=-i\ln{Z[J_{n}]}\;,$
$\displaystyle\Gamma[\langle\Phi_{n}\rangle_{J}]$ $\displaystyle=W[J_{n}]-\int
d^{4}x\,J_{n}(x)\langle\Phi_{n}\rangle_{J}(x)\;,$ (65)
$\displaystyle\mathcal{W}[J_{n},\phi_{n}]$
$\displaystyle=-i\ln{\mathcal{Z}[J_{n},\phi_{n}]}\;,$ $\displaystyle
H[\langle\phi^{\prime}_{n}\rangle_{J},\phi_{n}]$
$\displaystyle=\mathcal{W}[J_{n}]-\int
d^{4}x\,J_{n}(x)\langle\phi^{\prime}_{n}\rangle_{J}(x)\;,$ (66)

where $\langle\Phi_{n}\rangle_{J}$ and $\langle\phi^{\prime}_{n}\rangle_{J}$
represent the VEVs of $\Phi_{n}$ and $\phi^{\prime}_{n}$ in the presence of
sources $J_{n}$:

$\displaystyle\langle\Phi_{n}\rangle_{J}$ $\displaystyle\equiv\frac{\delta
Z}{\delta J_{n}}\;,$ $\displaystyle\langle\phi^{\prime}_{n}\rangle_{J}$
$\displaystyle\equiv\frac{\delta\mathcal{Z}}{\delta J_{n}}\;.$ (67)

Furthermore, by shifting the integration variable in $\mathcal{Z}$ it is
possible to conclude that these objects, although defined for different
generating functionals, are related to each other:

$\displaystyle\langle\phi^{\prime}_{n}\rangle_{J}$
$\displaystyle=\langle\Phi_{n}\rangle_{J}-\phi_{n}$ $\displaystyle\implies$
$\displaystyle H[\langle\phi^{\prime}_{n}\rangle_{J},\phi_{n}]$
$\displaystyle=\Gamma[\langle\phi^{\prime}_{n}\rangle_{J}+\phi_{n}]\;,$ (68)
$\displaystyle\phi_{n}$ $\displaystyle=\langle\Phi_{n}\rangle_{J}$
$\displaystyle\implies$ $\displaystyle H[0,\phi_{n}]$
$\displaystyle=\Gamma[\phi_{n}]\;.$ (69)

The last equation states that the quantum effective action of a theory can be
computed from its corresponding background effective action by turning off the
VEV of the fluctuations. Moreover, since one–particle irreducible (1PI)
diagrams are obtained by functional differentiation of the quantum effective
action one can conclude that: 1) $H[0,\phi_{n}]$ diagrams contain no external
lines of quantum fluctuations, 2) matrix elements in the original theory
(represented by $Z[J_{n}]$) can be computed by functionally differentiating
vacuum–to–vacuum diagrams in the background theory with respect to the
background fields, which is the result we required at the start of this
section. This last result is what justifies the fact that a product of
background fields, say $\overline{\psi}\sigma_{\mu\nu}\psi$, is related to the
composite operator $\overline{\Psi}\sigma_{\mu\nu}\Psi$.

There is an important caveat however: the gauge fixing procedure requires one
to be careful not to break the symmetries of the “original” $Z[J_{n}]$ when
building $\mathcal{Z}[J_{n},\phi_{n}]$ in order for these conclusions to hold.
In particular, we want to fix the gauge redundancy of the fluctuation:

$\delta^{\prime}A^{{}^{\prime}a}_{\mu}=\partial_{\mu}\epsilon^{a}(x)+g_{S}f^{abc}\epsilon^{c}(a_{\mu}^{b}+A^{{}^{\prime}b}_{\mu})\;,$
(70)

with matter fields transforming as usual and $\delta^{\prime}a_{\mu}^{a}=0$,
while keeping the lagrangian’s invariance under the background gauge
transformations:

$\delta a_{\mu}^{a}=\langle
0|\partial_{\mu}\epsilon^{a}+g_{S}f^{abc}A^{b}_{\mu}\epsilon^{c}|0\rangle=\partial_{\mu}\epsilon^{a}+g_{S}f^{abc}a_{\mu}^{b}\epsilon^{c}\;,\hskip
17.07182pt\delta A^{{}^{\prime}a}_{\mu}=\delta A^{a}_{\mu}-\delta
a_{\mu}^{a}=g_{S}f^{abc}A_{\mu}^{{}^{\prime}b}\epsilon^{c}\;.$ (71)

There is in fact only one gauge fixing term for the background generating
functional that fulfills such a requirement:
$-\frac{1}{2\xi}(\overline{D}^{\mu}A^{{}^{\prime}a}_{\mu})^{2}$, where
$\overline{D}^{\mu}A^{{}^{\prime}a}_{\mu}\equiv\partial^{\mu}A^{{}^{\prime}a}_{\mu}+g_{S}f^{abc}a^{b}_{\mu}A^{{}^{\prime}c\mu}$
acts as a background covariant derivative and $\xi$ is the arbitrary
gauge–fixing parameter.

Now that we have presented the theoretical framework of the separation of
fields that we will use, we are ready to obtain the computational tools that
we will need to build the OPE. Let us read in detail the Feynman rules from
the lagrangian $\mathcal{L}^{\text{HLbL}}$ after the separation of the fields
in classical background and quantum parts:

$\begin{split}\mathcal{L}^{\text{HLbL}}&=-\frac{1}{4}(f_{\mu\nu}^{a}+\overline{D}_{\mu}A^{{}^{\prime}a}_{\nu}-\overline{D}_{\nu}A^{{}^{\prime}a}_{\mu}+g_{S}f^{abc}A^{{}^{\prime}b}_{\mu}A^{{}^{\prime}c}_{\nu})^{2}-\frac{1}{4}(f_{\mu\nu}+\partial_{\mu}A^{{}^{\prime}a}_{\nu}-\partial_{\nu}A^{{}^{\prime}}_{\mu})^{2}\\\
&+(\overline{\psi}_{l}+\overline{\psi}^{\prime}_{l})(\\{i\not{\partial}-m\\}\delta_{lk}+e\hat{Q}\delta_{lk}\not{a}+g_{S}t^{a}_{lk}\not{a}^{a})(\psi_{l}+\psi^{\prime}_{l})\\\
&+(\overline{\psi}_{l}+\overline{\psi}^{\prime}_{l})(e\hat{Q}\delta_{lk}\not{A}^{\prime}+g_{S}t^{a}_{lk}\not{A}^{{}^{\prime}a})(\psi_{k}+\psi^{\prime}_{k})-\frac{1}{2\xi}\overline{D}^{\mu}A^{{}^{\prime}a}_{\mu}\overline{D}^{\mu}A^{{}^{\prime}a}_{\mu}-\mathcal{L}_{\text{ghosts}}\;,\end{split}$
(72)

where we have split the quark fields following the same (un)primed convention
and we have included colour indices $l$ and $k$. Besides, $f^{a}_{\mu\nu}$ and
$f_{\mu\nu}$ are background counterparts of the usual field strength tensors.

The previous lagrangian can be simplified very much. First, the parts with
only background field variables are constant from the point of view of the
path integral and therefore they may be adsorbed by its normalization
constant, as is usual for amplitudes of disconnected diagrams. In addition,
since HLbL scattering only involves vertices from strong interactions, the
quantum fluctuation part of the photon field is not relevant. In second place,
it is necessary to recall that the classical background field minimize the
action242424The classical background fields actually minimize the _quantum_
effective action, which in absence of external sources and to leading order is
equivalent to the same statement on the classical action. For quark and gluon
fields the vacuum expectation values of course do not receive perturbative
contributions since tadpole perturbative diagrams are zero at all orders.,
therefore the first functional derivative of the action with respect to each
field fluctuation vanishes, which means that terms of the lagrangian that are
linear in the field fluctuations do not contribute to the dynamics. Finally
for the gauge fixing parameter we choose $\xi=1$. In summary, the parts of the
lagrangian that are relevant to our case is:

$\begin{split}\mathcal{L}^{\text{HLbL}}&=-\frac{g_{S}}{2}f^{a\mu\nu}f^{abc}A^{{}^{\prime}b}_{\mu}A^{{}^{\prime}c}_{\nu}-\frac{1}{2}\Big{(}\overline{D}^{\mu}A^{{}^{\prime}a\nu}\overline{D}_{\mu}A^{{}^{\prime}a}_{\nu}+\overline{D}^{\mu}A^{{}^{\prime}a}_{\mu}\overline{D}^{\nu}A^{{}^{\prime}a}_{\nu}-\overline{D}^{\mu}A^{{}^{\prime}a\nu}\overline{D}_{\nu}A^{{}^{\prime}a}_{\mu}\Big{)}\\\
&-g_{S}f^{a\overline{b}\overline{c}}A^{{}^{\prime}\overline{b}\mu}A^{{}^{\prime}\overline{c}\nu}\overline{D}_{\mu}A^{{}^{\prime}a}_{\nu}-\frac{1}{4}g_{S}^{2}f^{abc}f^{a\overline{b}\overline{c}}A^{{}^{\prime}b}_{\mu}A^{{}^{\prime}c}_{\nu}A^{{}^{\prime}\overline{b}\mu}A^{{}^{\prime}\overline{c}\nu}\\\
&+\overline{\psi}^{\prime}_{l}(\\{i\not{\partial}-m\\}\delta_{lk}+e\hat{Q}\delta_{lk}\not{a}+g_{S}t^{a}_{lk}\not{a}^{a})\psi^{\prime}_{k}+\overline{\psi}_{l}(e\hat{Q}\delta_{lk}\not{A}^{\prime}+g_{S}t^{a}_{lk}\not{A}^{{}^{\prime}a})\psi^{\prime}_{k}\\\
&+\overline{\psi}^{\prime}_{l}(e\hat{Q}\delta_{lk}\not{A}^{\prime}+g_{S}t^{a}_{lk}\not{A}^{{}^{\prime}a})\psi_{k}+\overline{\psi}^{\prime}_{l}(e\hat{Q}\delta_{lk}\not{A}^{\prime}+g_{S}t^{a}_{lk}\not{A}^{{}^{\prime}a})\psi^{\prime}_{k}\\\
&-w^{*}_{a}\overline{D}_{\mu}\\{\overline{D}^{\mu}w_{a}+f_{abc}w_{c}A^{{}^{\prime}\mu}_{b}\\}\;,\end{split}$
(73)

where ghost fields are represented by $w$ and $w^{*}$ and the last line
corresponds to $\mathcal{L}_{\text{ghosts}}$, whose invariance under
background gauge transformations is guaranteed by the fact that ghost fields
transform like matter. As usual, the kernel of quadratic terms in a specific
path–integral variable is the inverse of corresponding free propagator, while
the rest are interaction vertices. These are very similar to the ones in a
theory with no background fields, the only difference being that any one line
can now be created or annihilated by the vacuum. Feynman rules for quark–gluon
vertices are summarized in figure 4. The vertices for gluon self interactions
can be found straightforwardly, but we do not quote them here because they are
not relevant at the perturbative level we have working at.

psiPsiPrimeAaPrime (5,5) i1,i2 $l$i2 o $a$o fermioni1,v,i2 gluonv,o arrow_len3mm decor.shape=circle,decor.size=1mmi1 |  | psiPrimePsiAaPrime (5,5) i1,i2 $k$i1 o $a$o fermioni1,v,i2 gluonv,o arrow_len3mm decor.shape=circle,decor.size=1mmi2 |
---|---|---|---
= $+ig_{S}\gamma^{\mu}t^{a}_{lk}\psi_{k}$ |  = $+ig_{S}\overline{\psi}_{l}t^{a}_{lk}\gamma^{\mu}$
| psiPrimePsiPrimeAaPrime (5,5) i1,i2 $l$i2 $k$i1 o $a$o fermioni1,v,i2 gluonv,o arrow_len3mm |  |
| = $+ig_{S}t^{a}_{lk}\gamma^{\mu}$ |
Figure 4: These figures show the three different types of quark interactions
with a quantum gluon. In the two diagrams at the top one quark line is
non–pertubatively annihilated by the vacuum. $l$ and $k$ represent the colour
of the quarks, $a$ represents the colour of the gluon and trivial quark flavor
indices are suppressed.

The quark free propagator is modified by the background fields in the usual
way:

$\displaystyle(\\{i\not{\partial}-m\\}\delta_{l^{\prime}l}\delta^{f^{\prime}f}+e\hat{Q}^{f^{\prime}f}\delta_{l^{\prime}l}\not{a}+g_{S}\delta^{f^{\prime}f}t^{a}_{l^{\prime}l}\not{a}^{a})S_{lk}^{fs}(x,y)$
$\displaystyle=i\delta^{4}(x-y)\delta^{f^{\prime}s}\delta_{l^{\prime}k}\;,$
(74)

therefore it can be computed recursively by supposing that the strength of the
background gauge fields is much smaller than the characteristic momentum of
the process of interest, which is a reasonable hypothesis in our context. Note
that the new indices $f^{\prime}f$ represent quark flavor and had been
suppressed previously. In the approximation of weak external electromagnetic
field, an expansion to $O(e)$ is enough. For the background gluon field it is
necessary to go to $O(g_{S}^{2})$ in agreement with the highest mass dimension
of the local operators that were chosen for the OPE.

In addition, we will see that due to the presence of a background the
propagator is not traslationally invariant. This of course also breaks
momentum conservation along the propagator in its Fourier transformed version.
Therefore there are in general three different types of momentum–space quark
propagators:

$\displaystyle S_{lk}^{fs}(p_{1},p_{2})$ $\displaystyle\equiv\int d^{4}x\int
d^{4}y\;e^{ip_{1}x}e^{-ip_{2}y}S_{lk}^{fs}(x,y)\;,$ (75) $\displaystyle
S_{lk}^{fs}(p_{1})$ $\displaystyle\equiv\int
d^{4}x\;e^{ip_{1}x}S_{lk}^{fs}(x,0)=\int\frac{d^{4}p_{2}}{(2\pi)^{4}}S_{lk}^{fs}(p_{1},p_{2})\;,$
(76) $\displaystyle\tilde{S}_{lk}^{fs}(p_{2})$ $\displaystyle\equiv\int
d^{4}x\;e^{-ip_{2}x}S_{lk}^{fs}(0,x)=\int\frac{d^{4}p_{1}}{(2\pi)^{4}}S_{lk}^{fs}(p_{1},p_{2})\;,$
(77)

By Taylor expanding the gauge fields we have intrinsically assumed that they
are soft, but not even that restores translation invariance of the propagator.
In summary, the free quark propagator in the momentum space (and in the
presence of background gluon and photon fields) at $O(eg^{2}_{S})$ is:

$\begin{split}S_{lk}^{fs}(p_{1},p_{2})&=(2\pi)^{4}\delta^{4}(p_{1}-p_{2})S^{0}_{p_{1}}\delta^{fs}\delta_{lk}\\\
&+iS^{0}_{p_{1}}\int_{q_{1}}\\{e\hat{Q}^{fs}\delta_{lk}\not{a}_{q_{1}}+g_{S}t^{a}_{lk}\delta^{fs}\not{a}^{a}_{q_{1}}\\}S^{0}_{p_{1}+q_{1}}(2\pi)^{4}\delta^{4}(p_{1}-p_{2}+q_{1})\\\
&-e\hat{Q}^{fs}g_{S}t^{a}_{lk}S^{0}_{p_{1}}\int_{q_{1}}\not{a}_{q_{1}}S^{0}_{p_{1}+q_{1}}\int_{q_{2}}\not{a}^{a}_{q_{2}}S^{0}_{p_{1}+q_{1}+q_{2}}(2\pi)^{4}\delta^{4}(p_{1}-p_{2}+q_{1}+q_{2})\\\
&-e\hat{Q}^{fs}g_{S}t^{a}_{lk}S^{0}_{p_{1}}\int_{q_{1}}\not{a}^{a}_{q_{1}}S^{0}_{p_{1}+q_{1}}\int_{q_{2}}\not{a}_{q_{2}}S^{0}_{p_{1}+q_{1}+q_{2}}(2\pi)^{4}\delta^{4}(p_{1}-p_{2}+q_{1}+q_{2})\\\
&-g_{S}^{2}t^{a}_{ll^{\prime}}t^{b}_{l^{\prime}k}\delta^{fs}S^{0}_{p_{1}}\int_{q_{1}}\not{a}^{a}_{q_{1}}S^{0}_{p_{1}+q_{1}}\int_{q_{2}}\not{a}^{b}_{q_{2}}S^{0}_{p_{1}+q_{1}+q_{2}}(2\pi)^{4}\delta^{4}(p_{1}-p_{2}+q_{1}+q_{2})\\\
&-ie\hat{Q}^{fs}g^{2}_{S}t^{a}_{ll^{\prime}}t^{b}_{l^{\prime}k}\Big{(}S^{0}_{p_{1}}\int_{q_{1}}\not{a}_{q_{1}}S^{0}_{p_{1}+q_{1}}\int_{q_{2}}\not{a}^{a}_{q_{2}}S^{0}_{p_{1}+q_{1}+q_{2}}\int_{q_{3}}\not{a}^{b}_{q_{3}}S^{0}_{p_{1}+q_{1}+q_{2}+q_{3}}\\\
&\hskip
68.28644pt+S^{0}_{p_{1}}\int_{q_{1}}\not{a}^{a}_{q_{1}}S^{0}_{p_{1}+q_{1}}\int_{q_{2}}\not{a}_{q_{2}}S^{0}_{p_{1}+q_{1}+q_{2}}\int_{q_{3}}\not{a}^{b}_{q_{3}}S^{0}_{p_{1}+q_{1}+q_{2}+q_{3}}\\\
&\hskip
68.28644pt+S^{0}_{p_{1}}\int_{q_{1}}\not{a}_{q_{1}}^{a}S^{0}_{p_{1}+q_{1}}\int_{q_{2}}\not{a}^{b}_{q_{2}}S^{0}_{p_{1}+q_{1}+q_{2}}\int_{q_{3}}\not{a}_{q_{3}}S^{0}_{p_{1}+q_{1}+q_{2}+q_{3}}\Big{)}(2\pi)^{4}\delta^{4}(p_{1}-p_{2}+\sum_{i=1,2,3}q_{i})\;.\end{split}$
(78)

Note that the Dirac delta is always under the effect of the integrals to its
left. Additionally, we use the convention:

$S^{0}(p)=i\frac{\not{p}+m}{p^{2}-m^{2}+i\epsilon}\equiv S^{0}_{p}\;.$ (79)

Note that we have kept the momentum–conservation delta of the $y$ vertex
explicitly for $S(p_{1},p_{2})$ in order to make the relation with $S(p)$ and
$\tilde{S}(p)$ more evident. At this point, we can use the expansions of the
gauge fields in (63). Since the local operators considered for the OPE contain
no derivatives of the photon fields and contain only up to one derivative of
the gluon field, it is enough to retain the first term of the expansion for
the photon field and the first two terms for the gluon. Therefore in the
momentum representation we can use the replacement252525Note that these
expressions were derived for soft insertions where momentum is leaving the
diagram. From the distributional point of view the derivatives are supposed to
act on test-functions [122].:

$\begin{split}a^{\mu}(q)&=\frac{i}{2}(2\pi)^{4}f^{\nu\mu}(0)\frac{\partial}{\partial
q^{\nu}}\delta^{4}(q)=-\frac{i}{2}(2\pi)^{4}\delta^{4}(q)f^{\nu\mu}(0)\frac{\partial}{\partial
q^{\nu}}\\\ a^{a\mu}(q)&=\frac{i}{2}f^{a\nu\mu}(0)\frac{\partial}{\partial
q^{\nu}}(2\pi)^{4}\delta^{4}(q)-\frac{1}{3}D^{\tau}f^{a\nu\mu}(0)\frac{\partial}{\partial
q^{\nu}}\frac{\partial}{\partial q^{\tau}}(2\pi)^{4}\delta^{4}(q)\\\
&=(2\pi)^{4}\delta^{4}(q)\Big{(}-\frac{i}{2}f^{a\nu\mu}(0)\frac{\partial}{\partial
q^{\nu}}-\frac{1}{3}D^{\tau}f^{a\nu\mu}(0)\frac{\partial}{\partial
q^{\nu}}\frac{\partial}{\partial q^{\tau}}\Big{)}\;,\end{split}$ (80)

which yields the following result:

$\begin{split}S_{lk}^{fs}(p_{1},p_{2})&=(2\pi)^{4}\Bigg{(}S^{0}_{p_{1}}\delta^{fs}\delta_{lk}\delta^{4}(p_{1}-p_{2})\\\
&+\frac{1}{2}\\{e\hat{Q}^{fs}\delta_{lk}f^{\mu_{1}\nu_{1}}(0)+g_{S}t^{a}_{lk}\delta^{fs}f^{a\mu_{1}\nu_{1}}(0)\\}\frac{\partial}{\partial
q^{\mu_{1}}_{1}}S^{0}_{p_{1}}\gamma_{\nu_{1}}S^{0}_{p_{1}+q_{1}}\delta^{4}(p_{1}-p_{2}+q_{1})\\\
&-\frac{i}{3}g_{S}t^{a}_{lk}\delta^{fs}\overline{D}^{\tau}f^{a\mu_{1}\nu_{1}}\frac{\partial}{\partial
q^{\tau}_{1}}\frac{\partial}{\partial
q^{\mu_{1}}_{1}}S^{0}_{p_{1}+q_{1}}S^{0}_{p_{1}}\gamma_{\nu_{1}}\delta^{4}(p_{1}-p_{2}+q_{1})\\\
&+\frac{1}{4}e\hat{Q}^{fs}g_{S}t^{a}_{lk}\Big{(}f^{\mu_{1}\nu_{1}}f^{a\mu_{2}\nu_{2}}+f^{a\mu_{1}\nu_{1}}f^{\mu_{2}\nu_{2}}\Big{)}\\\
&\hskip 28.45274pt\times\frac{\partial}{\partial
q^{\mu_{1}}_{1}}\frac{\partial}{\partial
q^{\mu_{2}}_{2}}\Big{(}S^{0}_{p_{1}}\gamma_{\nu_{1}}S^{0}_{p_{1}+q_{1}}\gamma_{\nu_{2}}S^{0}_{p_{1}+q_{1}+q_{2}}\delta^{4}(p_{1}-p_{2}+q_{1}+q_{2})\Big{)}\Bigg{|}_{q_{1,2}=0}\\\
&+\frac{1}{4}g_{S}^{2}t^{a}_{ll^{\prime}}t^{b}_{l^{\prime}k}\delta^{fs}\Big{(}f^{a\mu_{1}\nu_{1}}f^{b\mu_{2}\nu_{2}}\Big{)}\\\
&\hskip 28.45274pt\times\frac{\partial}{\partial
q^{\mu_{1}}_{1}}\frac{\partial}{\partial
q^{\mu_{2}}_{2}}\Big{(}S^{0}_{p_{1}}\gamma_{\nu_{1}}S^{0}_{p_{1}+q_{1}}\gamma_{\nu_{2}}S^{0}_{p_{1}+q_{1}+q_{2}}\delta^{4}(p_{1}-p_{2}+q_{1}+q_{2})\Big{)}\Bigg{|}_{q_{1,2}=0}\\\
&+\frac{1}{8}e\hat{Q}^{fs}g^{2}_{S}t^{a}_{lk^{\prime}}t^{b}_{k^{\prime}k}\Big{(}f^{\mu_{1}\nu_{1}}f^{a\mu_{2}\nu_{2}}f^{b\mu_{3}\nu_{3}}+f^{a\mu_{1}\nu_{1}}f^{\mu_{2}\nu_{2}}f^{b\mu_{3}\nu_{3}}+f^{a\mu_{1}\nu_{1}}f^{b\mu_{2}\nu_{2}}f^{\mu_{3}\nu_{3}}\Big{)}\\\
&\hskip 28.45274pt\times\frac{\partial}{\partial
q^{\mu_{1}}_{1}}\frac{\partial}{\partial
q^{\mu_{2}}_{2}}\frac{\partial}{\partial
q^{\mu_{3}}_{3}}\Big{(}S^{0}_{p_{1}}\gamma_{\nu_{1}}S^{0}_{p_{1}+q_{1}}\gamma_{\nu_{2}}S^{0}_{p_{1}+q_{1}+q_{2}}\gamma_{\nu_{3}}S^{0}_{p_{1}+q_{1}+q_{2}+q_{3}}\delta^{4}(p_{1}-p_{2}+\sum_{i=1,2,3}q_{i})\Big{)}\Big{|}_{q_{1,2,3}=0}\Bigg{)}\;.\end{split}$
(81)

This expression for the free quark propagator can be understood as an
expansion in terms of diagrams (see figure 5) with increasing number of
(background) gauge bosons.

quarkProp (5,5) i o heavyi,o arrow_len3mm |  | quarkProp0 (5,5) i o fermioni,o arrow_len3mm |  | quarkPropE (5,5) i o t fermioni,v,o photont,v decor.shape=circle,decor.size=1mmt arrow_len3mm |  | quarkPropG (5,5) i o t fermioni,v,o gluont,v decor.shape=circle,decor.size=1mmt arrow_len3mm
---|---|---|---|---|---|---
= | + | +
|  | quarkPropEG (5,5) i o t1,t2 fermioni,v1,v2,o photont1,v1 gluont2,v2 decor.shape=circle,decor.size=1mmt1,t2 arrow_len3mm |  | quarkPropGE (5,5) i o t1,t2 fermioni,v1,v2,o gluon,wiggly_len=4mmt1,v1 photont2,v2 decor.shape=circle,decor.size=1mmt1,t2 arrow_len3mm |  | quarkPropGG (5,5) i o t1,t2 fermioni,v1,v2,o gluont1,v1 gluont2,v2 decor.shape=circle,decor.size=1mmt1,t2 arrow_len3mm
| + | + | +
|  | quarkPropEGG (5,5) i o t1,t2,t3 fermioni,v1,v2,v3,o photont1,v1 gluont2,v2 gluont3,v3 decor.shape=circle,decor.size=1mmt1,t2,t3 arrow_len3mm |  | quarkPropGEG (5,5) i o t1,t2,t3 fermioni,v1,v2,v3,o gluon,wiggly_len=4mmt1,v1 photont2,v2 gluont3,v3 decor.shape=circle,decor.size=1mmt1,t2,t3 arrow_len3mm |  | quarkPropGGE (5,5) i o t1,t2,t3 fermioni,v1,v2,v3,o gluon,wiggly_len=4mmt1,v1 gluont2,v2 photont3,v3 decor.shape=circle,decor.size=1mmt1,t2,t3 arrow_len3mm
| + | + | +
Figure 5: This figure shows the expansion of the free quark propagator in a
background of gauge fields in terms of diagrams with interactions with gluons
and photons that are created/annihilated in the vacuum. The order in which
diagrams appear in the sum corresponds to the order of terms in equation (81).

The “free” gluon propagator contains interactions with the background gluon
fields in a similar but more involved way than the quark one. However, at the
order that we are interested in there appear no such propagators, so we refer
the interested reader to see the details of the computation in section 2 of
[123].

##  4 Computation of un–renormalized Wilson coefficients

In the previous section we obtained expressions for the quark and gluon
fluctuations propagators which contained background insertions of vacuum
expectation values (VEV) such as $f^{a\mu\nu}$ and $f^{\mu\nu}$. Furthermore
we saw that vertices from the Dyson series also introduce VEVs of quark
operators, thus giving us all the tools required to build the OPE of
$\Pi^{\mu_{1}\mu_{2}\mu_{3}}$ with background fields and find the Wilson
coefficients that require to compute
$\Pi_{F}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\mu_{5}}$. Concerning the actual
computation of Wilson coefficients, let us start by considering the one
related to $S_{1,\mu\nu}=ee_{f}F_{\mu\nu}$. This term represents the
configuration in which hard momenta travels through all internal lines of the
diagrams, thus, there are no cut lines. The leading order contribution for
this configuration is given by the quark loop (see figure 6), where different
contributions are obtained by inserting the soft photon in different sides of
the triangle and/or inverting the orientation of the loop. Since
$S_{1,\mu\nu}$ is the operator with the lowest dimension in the OPE, its
Wilson coefficient is expected to give the most relevant contribution to
$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\mu_{5}}_{F}$ and therefore to $a_{\mu}$.
When splitting the quark fields in the currents of
$\Pi^{\mu_{1}\mu_{2}\mu_{3}}$, this diagram comes from the term that contains
only quantum fluctuations. In the end, the contribution from the Wilson
coefficient of $S_{1,\mu\nu}$ to
$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\mu_{5}}_{F}$:

$\begin{split}\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\nu_{4}}_{F(S_{1})}&=i\frac{N_{c}}{2}\int\frac{d^{4}p}{(2\pi)^{4}}\sum_{f}e_{f}^{4}\frac{\partial}{\partial
q_{4\nu_{4}}}\sum_{\sigma(1,2,4)}\operatorname{Tr}\Big{\\{}\gamma^{\mu_{3}}S^{0}(p+q_{1}+q_{2}+q_{4})\gamma^{\mu_{4}}\\\
&\hskip 28.45274pt\times
S^{0}(p+q_{1}+q_{2})\gamma^{\mu_{1}}S^{0}(p+q_{2})\gamma^{\mu_{2}}S^{0}(p)\Big{\\}}\Bigg{|}_{q_{4}=0}\;,\end{split}$
(82)

where $N_{c}$ is the number of quark colors, $f$ represents quark flavor and
$\sigma(1,2,4)$ represents a permutation over the set
$\\{(q_{i},\mu_{i})|\,i\in\\{1,2,4\\}\\}$. The derivative with respect to the
soft photon momentum can be traced back to (81). Note that propagator depends
implicitly on the quark flavor through the masses, assumed to be all equal for
the light flavours considered. Besides, the effect of the derivative on the
propagators is to duplicate them:

$\lim_{q_{4}\rightarrow 0}\frac{\partial}{\partial
q_{\nu_{4}}}S(p+q_{4})=i\lim_{q_{4}\rightarrow
0}S^{0}(p+q_{4})\gamma^{\nu_{4}}S(p+q_{4})=iS^{0}(p)\gamma^{\nu_{4}}S^{0}(p)\;.$
(83)

quarkLoop (5,5) i1,i2 o1,o2,o3 phantomv1,v2,v3,v4,v5 phantomv5,v6,v7,v8,v9
phantomv9,v10,v11,v12,v1 photoni1,v1 photonv9,i2 photonv5,o2 photonv3,o1
plainv1,v2,v3,v4,v5 fermionv5,v9,v1 arrow_len3mm
decor.shape=circle,decor.size=1mmo1 (-0.3w,0.2h)o1

Figure 6: Representative diagram of the leading order contribution to the
Wilson coefficient of $S_{1,\mu\nu}$ in the OPE of
$\Pi^{\mu_{1}\mu_{2}\mu_{3}}$. The black dot represents creation/annihilation
of a line by the background fields in the vacuum.

The focus of our work is on the contribution from $S_{1,\mu\nu}$, which gives
a much larger contribution to $a_{\mu}$ than any other operator, so we will
not discuss the contributions from other operators in detail. Instead, we give
a brief overview of the computation. We refer the interested reader to [56].

Contributions with one cut quark line and at most one soft gauge boson
insertion ($S_{2-5,\mu\nu}$ and $S_{7,\mu\nu}$) are obtained at leading order
from the diagrams in figure 7. Their corresponding amplitudes are computed
from terms in $\Pi^{\mu_{1}\mu_{2}\mu_{3}}$ that contain two soft quark fields
and require no vertices from the Dyson series expansion. Soft gluon or photon
insertions on quark hard lines, if necessary, come from propagators of quark
fluctuations as seen in the previous section.

cutQuarkLine (5,5) i1,i2 o1,o2,o3 phantomv1,v2,v3,v4,v5 phantomv5,v6,v7,v8,v9 phantomv9,v10,v11,v12,v1 photoni1,v1 photonv9,i2 photonv5,o2 photonv11,o1 fermionv1,v5,v9 plainv9,v10 plainv12,v1 arrow_len3mm decor.shape=circle,decor.size=1mmo1,v10,v11,v12 (-0.9w,0.5h)o1 (0,-0.06h)v10 (0,0.06h)v12 |  cutQuarkLineGluon (5,5) i1,i2 o1,o2,o3 phantomv1,v2,v3,v4,v5 phantomv5,v6,v7,v8,v9 phantomv9,v10,v11,v12,v1 photoni1,v1 photonv9,i2 photonv5,o2 photonv11,o3 gluonv3,o1 plainv1,v2,v3,v4,v5 fermionv5,v9 plainv9,v10 plainv12,v1 arrow_len3mm decor.shape=circle,decor.size=1mmo1,v10,v11,v12,o3 (-0.3w,0.2h)o1 (-0.9w,-0.5h)o3 (0,-0.06h)v10 (0,0.06h)v12 |  cutQuarkLinePhoton (5,5) i1,i2 o1,o2,o3 phantomv1,v2,v3,v4,v5 phantomv5,v6,v7,v8,v9 phantomv9,v10,v11,v12,v1 photoni1,v1 photonv9,i2 photonv5,o2 photonv3,o1 plainv1,v2,v3,v4,v5 fermionv5,v9 plainv9,v10 plainv12,v1 arrow_len3mm decor.shape=circle,decor.size=1mmo1,v10,v12 (-0.3w,0.2h)o1 (0,-0.075h)v10 (0,0.075h)v12
---|---|---
Figure 7: Representative diagrams of the leading order contribution to the
Wilson coefficient of $S_{2,\mu\nu}$ (first diagram), $S_{3,4,7,\mu\nu}$
(second diagram) and $S_{5,\mu\nu}$ (third diagram) in the OPE of
$\Pi^{\mu_{1}\mu_{2}\mu_{3}}$. The black dot represents creation/annihilation
of a line by the background fields in the vacuum.

Let us now consider the operator with two cut gluon lines, that is,
$S_{6,\mu\nu}$. Diagrams contributing to this operator are very similar to the
quark loop of $S_{1,\mu\nu}$, but they have two soft gluon insertions (see
figure 8). As with the first quark loop, these insertions must be permuted in
all possible ways to obtain the full contribution.

twoCutGluonQuarkLoop (5,5) i1,i2,i3 o1,o2,o3 phantomv1,v2,v3,v4,v5
phantomv5,v6,v7,v8,v9 phantomv9,v10,v11,v12,v1 photoni1,v1 photonv9,i3
photonv5,o2 gluonv3,o1 gluonv7,o3 photonv11,i2 plainv1,v3 fermionv3,v5
plainv5,v7 fermionv7,v9 plainv9,v11 fermionv11,v1 arrow_len3mm
decor.shape=circle,decor.size=1mmo1,o3,i2 (0.025w,0h)i2 (-0.3w,-0.2h)o3
(-0.3w,0.2h)o1

Figure 8: Representative diagram of the leading order contribution to the
Wilson coefficient of $S_{6,\mu\nu}$ in the OPE of
$\Pi^{\mu_{1}\mu_{2}\mu_{3}}$. The black dot represents creation/annihilation
of a line by the background fields in the vacuum.

Let us now consider operators with four quark background insertions
($S_{8,\mu\nu}$). Diagrams that contribute to the Wilson coefficients of this
operator correspond to the quark loop with two cut quark lines, therefore the
diagram is divided in two parts, which have to be connected by a gluon (see
figure 9). There are six different ways in which the virtual gluon line can
connect the two parts of the diagram and all have to be accounted for. The
corresponding two gluon–quark vertices are responsible for the $\alpha_{S}$
coefficient of $S_{8,\mu\nu}$.

twoCutQuarkLines (5,5) i1,i2 o1,o2,o3 phantomv1,v2,v3,v4,v5
phantomv5,v6,v7,v8,v9 phantomv9,v10,v11,v12,v1 photoni1,v1 photonv9,i2
photonv5,o2 gluon,wiggly_len=7mm,right=0.2,tension=0.8v3,v11 photonv7,o3
fermionv1,v3 fermionv3,v5 plainv5,v6 plainv8,v9,v10 plainv12,v1 arrow_len3mm
decor.shape=circle,decor.size=1mmv6,v7,v8,v10,v12,o3 (0,-0.075h)v10
(0,0.075h)v12 (-0.3w,-0.15h)o3 (0.060w,-0.040h)v8 (-0.060w,0.030h)v6
(0,0.160h)v11 curly_len2mm
---
Figure 9: Representative diagram of the leading order contribution to the
Wilson coefficient of $S_{8,\mu\nu}$ in the OPE of
$\Pi^{\mu_{1}\mu_{2}\mu_{3}}$. The black dot represents creation/annihilation
of a line by the background fields in the vacuum.

Except for $S_{1,\mu\nu}$, contributions from all operators to
$\Pi_{F}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\mu_{5}}$ depend on the susceptibilities
$X_{i}^{S}$. By definition these are non–perturbative quantities which are
usually computed either by lattice, models and/or educated guesses. The most
well–known one is $X_{5}$,262626Note that we have suppressed the $S$ index.
This was done because, as we will see in the next section, the elements of the
OPE need renormalization and therefore a new set of susceptibilities $X_{i}$
is defined in terms of the renormalized operators. because it is related to
the quark condensate which is a common subject of study in lattice
computations. A more recent version of the review cited in [56] can be found
in [124], where figure 14, table 22 and references therein represent a
thorough compilation of results for the quark condensate. The rest of the
susceptibilities are not so well–known and their numerical values are
estimated in [56] by a combination of models and educated guesses.

Up to this point we have presented all unrenormalized Wilson coefficients
associated with operators in (59) and, more importantly, their contribution to
$\partial^{\mu_{5}}\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$. Computation of the
Wilson coefficients is however not yet complete, for renormalization of the
OPE elements has not been taken into account. In contrast to the usual
situation in perturbative computations, we have not encountered ultraviolet
divergences in the Wilson coefficients of this section. In fact, except for
$S_{1,\mu\nu}$ and $S_{6,\mu\nu}$ all of their leading order contributions are
at tree level. As we will see in the next sections for the quark loop, the
Wilson coefficients of these two operators, although finite, have infrared
contributions that are renormalized by the quark masses. Such singularities
scale as logarithms and negative powers of $m_{f}$. These singular terms are
problematic in a twofold way. From a computational perspective these singular
factors may spoil convergence of the perturbative computation when the momenta
of the process, namely $Q_{i}$, get much bigger that the mass scale of the
quarks, which is actually our situation. From a conceptual point of view it is
also questionable to have Wilson coefficients with infrared contributions: in
the OPE framework they are meant to represent the contribution from the parts
of the diagram through which the external very high momenta travel. In the
next section we will present how renormalization of the the product of
background fields “cures” these infrared divergences and thus completes the
separation of low and high energy contributions of the OPE.

##  5 OPE of $\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ in an electromagnetic
background field: Renormalization

In this section we will present the renormalization program for the operators
that form the OPE for $\Pi^{\mu_{1}\mu_{2}\mu_{3}}$ in the $\overline{MS}$
scheme.

The Wilson coefficients that were presented in the previous section are UV
finite at the computed order, but they do have infrared divergent terms such
as $1/m_{f}^{2}$ and $\ln\\{Q^{2}_{i}/m_{f}^{2}\\}$ that are regularized by
the quarks masses. These terms may spoil the convergence of the perturbative
expansion. Moreover, the Wilson coefficients should not have infrared
contributions in the first place, therefore it should be possible to safely to
compute them in the massless quarks limit. There is an additional kind of low
energy effects that may affect Wilson coefficients: the ones arising from
diagrams where soft quark and gluon272727This does not apply for photon soft