LU-TP 22-63

November 2022

Constraints on the hadronic light-by-light in the Melnikov-Vainshtein regime
Johan Bijnens ${}^{a,\,\dagger}$, Nils Hermansson-Truedsson
${}^{a,\,b,\,\ddagger}$, Antonio Rodríguez-Sánchez ${}^{c,\,\ast}$
a Department of Astronomy and Theoretical Physics, Lund University,

Box 43, SE 221 00 Lund, Sweden

b Albert Einstein Center for Fundamental Physics, Institute for Theoretical
Physics, Universität Bern, Sidlerstrasse 5, 3012 Bern, Switzerland

c Scuola Internazionale Superiore di Studi Avanzati (SISSA),

via Bonomea, 265 - 34136 Trieste, Italy

The muon anomalous magnetic moment continues to attract attention due to the
possible tension between the experimentally measured value and the theoretical
Standard Model prediction. With the aim to reduce the uncertainty on the
hadronic light-by-light contribution to the magnetic moment, we derive short-
distance constraints in the Melnikov-Vainshtein regime which are useful for
data-driven determinations. In this kinematical region, two of the four
electromagnetic currents are close in the four-point function defining the
hadronic light-by-light tensor. To obtain the constraints, we develop a
systematic operator product expansion of the tensor in question to next-to-
leading order in the expansion in operators. We evaluate the leading in
$\alpha_{s}$ contributions and derive constraints for the next-to-leading
operators that are also valid nonperturbatively.

†<EMAIL_ADDRESS>
‡<EMAIL_ADDRESS>
∗<EMAIL_ADDRESS>
###### Contents

  1. 1 Introduction
  2. 2 Generalities of the hadronic light-by-light
  3. 3 The two-current OPE
  4. 4 Perturbative regime
    1. 4.1 Dimension $D=3$
    2. 4.2 Dimension $D=4$
      1. 4.2.1 Dimension $D=4$ quark current operator
      2. 4.2.2 Dimension $D=4$ photon field strength operator
      3. 4.2.3 Expanded dimension $D=3$ operator
      4. 4.2.4 Final $\hat{\Pi}_{i}$ at dimension $D=4$
      5. 4.2.5 Renormalization of the OPE
  5. 5 Nonperturbative regime
    1. 5.1 Dimension $D=3$
    2. 5.2 Dimension $D=4$
      1. 5.2.1 Octet case
      2. 5.2.2 Singlet case
  6. 6 Conclusions

##  1 Introduction

The object of study in this work is the hadronic light-by-light (HLbL)
tensor111The integral over $q_{4}$ simply removes the delta function for
conservation of momenta.,

$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}=-i\int\frac{d^{4}q_{4}}{(2\pi)^{4}}\left(\prod_{i=1}^{4}\int
d^{4}x_{i}\,e^{-iq_{i}x_{i}}\right)\langle
0|T\left(\prod_{j=1}^{4}J^{\mu_{j}}(x_{j})\right)|0\rangle\,.$ (1)

Above, the electromagnetic current is defined as

$J^{\mu}(x)={\sum_{q}}e_{q}\,\bar{q}(x)\gamma^{\mu}q(x)\,,$ (2)

with $q=u,d,s$ and $e_{q}$ the associated light-quark charge matrix,
$e_{q}=\frac{1}{3}\left(2,-1,-1\right)$. In spite of its very simple
definition, only involving color symmetry and three light quarks, its emerging
dynamics is extremely complex and rich, as a consequence of the underlying
nonperturbative color dynamics. Its multi-scale nature represents a
theoretical challenge to understand where and how a picture in terms of
approximately free quarks and gluons can be used to describe it, and how to
separate the corresponding asymptotic QCD dynamics from the nonperturbative
hadronic one. Achieving such a task with a systematic and unambiguous power
counting, employing rigorous OPE methods, is one of the motivations of this
research line, started in Refs. [1, 2, 3] and continued here in a different
kinematic regime, all of them in the region where one of the external photons
becomes soft. The motivation to focus on that region is the key
phenomenological role played by the HLbL tensor in theoretically assessing the
anomalous magnetic moment of the muon via the diagram depicted in Fig. 1.
Indeed, its contribution represents one of the main uncertainties in the most
recent Standard Model (SM) value given in the White Paper [4],

$a_{\mu}^{\textrm{WP}}=\frac{(g-2)_{\mu}}{2}=116\,591\,810(43)\times
10^{-11}\,.$ (3)

The quoted value is $4.2\sigma$ away from the experimental one [5, 6],

$a_{\mu}^{\textrm{exp}}=116\,592\,061(41)\times 10^{-11}\,.$ (4)

However, a recent lattice evaluation [7] of the leading hadronic vacuum
polarization contribution (HVP) suggests that $R$-ratio based evaluations
might be underestimating its contributions, possibly explaining the tension
without any need of Beyond the Standard Model (BSM). A short discussion and
more relevant references can be found in Ref. [8].

While clarifying the situation with respect to the HVP side has become the
most urgent matter, further improvements on the HLbL evaluation are needed in
order to match the improving experimental precision. The first proper
evaluations date back to 1990s [9, 10, 11, 12] and were mainly based on
models. Our study is part of the data-driven evaluation of HLbL, which is
based on dispersion relations for the leading low-energy contributions [13,
14]. The subleading, but not negligible, contributions from the short- and
intermediate-distance regions, i.e., involving any large or intermediate loop
momenta, are quite challenging to assess within the dispersive approach, and
one needs to rely on models and other constraints. A thorough discussion can
be found in the White Paper [4]. Here we simply quote the dispersive White
Paper prediction [4]

$\displaystyle a_{\mu}^{\textrm{HLbL}}=92(19)\times 10^{-11}\,,$ (5)

which is based on Refs. [15, 16, 14, 17, 18, 1, 19, 20, 21, 22, 23, 24,
25].222Let us also point out that major improvements on the lattice-based
evaluation have been made since them. See Ref. [26].

Systematic uncertainties from the short-distance (SD) behaviour of the HLbL
tensor can be controlled with the help of short-distance constraints (SDCs).
Such nontrivial information about the analytical structure of the HLbL tensor
can be obtained through rigorous Operator Product Expansion (OPE) methods.
These OPEs are constructed for different limits of the Euclidean virtualities
$Q_{i}^{2}=-q_{i}^{2}$ associated to the virtual photons in Fig. 1. In Refs.
[1, 2] we showed how a systematic OPE [27, 28], different from the vacuum OPE,
which is valid when all Euclidean loop momenta of the HLbL tensor are large
[29], can be applied to the HLbL when all the loop momenta in Fig. 1 are
large, i.e. $Q_{i}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$, in spite of the soft
nature ($q_{4}\rightarrow 0$) of the photon coupling to the external magnetic
field. The naive massless quark loop happened to be the leading contribution,
directly related to the fact that it comes from the photon operator
$F_{\mu\nu}$, that, up to higher orders in QED, is not affected by quark-gluon
interactions. The perturbative corrections to the Wilson coefficient of this
operator were calculated in Ref. [3].

Figure 1:  The HLbL contribution to $a_{\mu}$. The grey blob contains hadronic
contributions. The photon with momentum $q_{4}$ is the external magnetic
field. The bottom line is the muon and the wiggly lines are photons.

Another regime of short-distances for the HLbL tensor and its contribution to
the muon $g-2$ is where two of the Euclidean momenta are large. The study of
this area was pioneered by Ref. [15]. In this limit a second Euclidean loop
momentum, apart from the static $q_{4}\rightarrow 0$ photon corresponding to
the magnetic field, becomes much smaller than the other two. The leading term
and some of the implications were worked out in Ref. [15]. While there is some
debate about the phenomenological consequences of the corresponding SD
constraints [30, 19, 31], their validity is well-established and relevant for
constraining the HLbL, see e.g. Refs. [31, 19, 30, 32, 33, 34, 35, 36, 37, 38,
39, 40, 41] for implementations of the SD constraints for various models and
intermediate contibutions. The White Paper estimate for the short-distance
contribution to the HLbL is [4]

$\displaystyle\Delta a_{\mu}^{\textrm{HLbL},\,\textrm{SDC}}=15(10)\times
10^{-11}\,.$ (6)

The main goal of this work consists in developing the associated OPE beyond
the limits obtained in [15], adding power corrections and parameterizing the
nonperturbative ignorance in a few form factors whose energy dependence are
known for large momenta.

In Sec. 2 we present some more details on the general properties of the HLbL.
The OPE studied in this work is formulated in Sec. 3 through next-to-leading
order (NLO). In Sec. 4 we focus on the region where the matrix elements can be
computed perturbatively up to NLO in the energy expansion and compare with the
results of Ref. [3]. The full agreement with the corresponding expansion of
the OPE of Ref. [3] is also shown in that section, leading to a very strong
consistency check of the interplay of the different two regions. Finally in
Sec. 5 we perform a nonperturbative extrapolation at leading order in
$\alpha_{s}$ at the hard scale, parametrizing the nonperturbative ignorance in
a few scalar functions depending on only a single squared momentum.
Conclusions are then given in Sec. 6.

##  2 Generalities of the hadronic light-by-light

The HLbL contribution to the muon magnetic moment, $a_{\mu}^{\textrm{HLbL}}$,
is depicted diagrammatically in Fig. 1. The underlying object to consider is
thus a four-point function of electromagnetic currents, i.e. the HLbL tensor
as defined in Eq. (1). We denote the virtual photon momenta by $q_{1}$,
$q_{2}$, $q_{3}$, and the external photon momentum $q_{4}$ which is soft for
the $(g-2)_{\mu}$ kinematics, it corresponds to the external magnetic field.
These momenta are chosen to satisfy $q_{1}+q_{2}+q_{3}+q_{4}=0$. We denote the
magnitudes of the three momenta in the Euclidean by $Q_{i}$, i.e.
$q_{i}^{2}=-Q_{i}^{2}$. The tensor in (1) satisfies the Ward identities

$\displaystyle q_{i,\,\mu_{i}}\,\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$
$\displaystyle=0\,,\qquad i=1,\dots,4\,,$ (7)

which immediately leads to

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}=-q_{4,\,\nu_{4}}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4,\,\mu_{4}}}\,.$ (8)

The above relation implies that the derivative of the HLbL tensor at $q_{4}=0$
can be used to fully determine the HLbL contribution to the magnetic moment,
$a_{\mu}^{\textrm{HLbL}}$ [42].

For this reason, we will in the following discuss the derivative of the HLbL
tensor

$\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4,\mu_{4}}}\,.$ (9)

This tensor can be decomposed into a set of $19$ independent scalar form
factors at $q_{4}\to 0$ [15, 13]. However, these may be generated from only
$5$ different form factors $\hat{\Pi}_{j}$ by using crossing relations [13].
The precise relation to the derivative can be found in Ref. [13].

As explained in detail in Refs. [13, 14], from the five independent
$\hat{\Pi}_{j}$ functions one can construct twelve scalar functions
$\overline{\Pi}_{i}$ such that $a_{\mu}^{\textrm{HLbL}}$ is given by

$\displaystyle a_{\mu}^{\textrm{HLbL}}=\frac{2\alpha^{3}}{3\pi^{2}}$
$\displaystyle\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}\sum_{i=1}^{12}T_{i}(Q_{1},Q_{2},\tau)\,\overline{\Pi}_{i}(Q_{1},Q_{2},\tau)\,,$
(10)

where the $T_{i}$ are known kernel functions and $\tau$ is defined by
$Q_{3}^{2}=Q_{1}^{2}+Q_{2}^{2}+2\tau Q_{1}Q_{2}$.

##  3 The two-current OPE

The HLbL tensor in Eq. (1) may be rewritten as

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}=$
$\displaystyle\sum_{j,k}\frac{ie_{q_{j}}e_{q_{k}}}{e^{2}}\int\frac{d^{4}q_{4}}{(2\pi)^{4}}\int
d^{4}x_{1}\int d^{4}x_{2}\,e^{-i(q_{1}x_{1}+q_{2}x_{2})}$
$\displaystyle\times\langle
0|T(J_{j}^{\mu_{1}}(x_{1})J_{k}^{\mu_{2}}(x_{2}))|\gamma^{\mu_{3}}(q_{3})\gamma^{\mu_{4}}(q_{4})\rangle\,,$
(11)

where we introduced a flavor-index on the current such that

$J^{\mu}_{j}=\bar{q}_{j}\gamma^{\mu}q_{j}\,.$ (12)

Note that two of the photons have been captured in the external state, which
can be contracted with electromagnetic interaction terms in the implicit Dyson
series with the QCD+QED action. The Lorentz indices on the photons in the
external state is a convenient short-hand notation for contraction with a
photon field $A^{\alpha}$ in the Dyson series according to

$A^{\alpha}(x)|\gamma^{\mu}(q)\rangle\equiv g^{\alpha\mu}e^{-iqx}\,.$ (13)

The difference to the usual contraction

$A^{\alpha}(x)|\gamma(q)\rangle=\epsilon^{\alpha}(q)e^{-iqx}\,,$ (14)

is that we factor out the photon polarization vector $\epsilon^{\alpha}(q)$.
The regime we study in this paper is where we can perform a proper OPE on the
time ordered product of two quark currents,
$T(J_{j}^{\mu_{1}}(x_{1})J_{k}^{\mu_{2}}(x_{2}))$. From the operators
apprearing in this OPE we only need to keep those that can mediate the
transition to two extra photons, one of which, $\gamma^{\mu_{4}}(q_{4})$, is
in the static $q_{4}\rightarrow 0$ limit relevant for the $(g-2)_{\mu}$.

The OPE is valid as far as there is an energy scale flowing between both
space-time points much larger than the rest of scales entering into the
process, which can then be factored out from the long-distance problem. Often
the mass of a heavy propagator plays such a role and we refer to the
corresponding OPE application as integrating out a heavy particle. In this
case the propagators, as far as the quark-gluon picture holds, are (almost)
massless. Nevertheless, in the kinematic regime in which most of the momentum
$q_{1}$ (which must be large enough with respect to the dynamically generated
$\Lambda_{\mathrm{QCD}}$ scale, in such a way that the quark-gluon picture
holds) that enters at $x_{1}$, leaves the process at $x_{2}$, an OPE can be
performed, taking $x_{1}$ and $x_{2}$ close to each other. We thus define
$\hat{q}$ via

$\displaystyle\hat{q}$ $\displaystyle=\frac{1}{2}\left(q_{1}-q_{2}\right)\,,$
$\displaystyle q_{1,2}=\pm\hat{q}-\frac{1}{2}\left(q_{3}+q_{4}\right)\,.$ (15)

The regime we consider is taking $q_{3}$ and $q_{4}$ much smaller than
$\hat{q}$. The large scale is then given by $\hat{Q}=\sqrt{-\hat{q}^{2}}$. Of
course, in the integration over loop momenta in Eq. (10) one has to consider
also the cases with $q_{1}$ or $q_{2}$ small, but these can be handled
analogously.

When dealing with the contribution to $a_{\mu}^{\textrm{HLbL}}$ more
convenient variables are $\overline{Q}_{3}=Q_{1}+Q_{2}$,
$\delta_{12}=Q_{1}-Q_{2}$ and $Q_{3}$. The sizes of these variables are
constrained, when $q_{4}\to 0$, by $q_{1}+q_{2}+q_{3}=0$ implying
$|\delta_{12}|\leq Q_{3}$. The large $\hat{q}$ and the large
$\overline{Q}_{3}$ expansions are directly related. We will use the $\hat{Q}$
large for the OPE and the large $\overline{Q}_{3}$ for the final expressions
for the $\hat{\Pi}_{i}$.

Since the object in Eq. (3) has dimension $-2$ in energy, a local operator of
dimension $D$ is going to give a contribution suppressed by
$\frac{1}{\overline{Q}_{3}^{D-2}}$. Let us now perform the derivation of the
OPE up to $D=4$. We may start by the contributions given by operators
involving quark fields. Let us first simplify the notation of the integration
and matrix element in Eq. (3). To do this, we note that at leading order in
$\alpha_{s}(\overline{Q}_{3})$ one has to contract two of the quark fields in
the two electromagnetic currents in the product
$J_{j}^{\mu_{1}}(x_{1})J_{k}^{\mu_{2}}(x_{2})$, thus giving rise to a quark
propagator and the restriction $j=k$. The only way of achieving a potential
$j\neq k$ contribution at NLO would be with a gluon propagator between the two
quark lines. In order to get a local contribution of $D\leq 4$ (terms of
increasing dimension in the OPE are suppressed by larger powers of
$\frac{Q_{3},\Lambda_{\mathrm{QCD}}}{\overline{Q}_{3}}$) one needs to close at
least one extra quark line, but then the associated color trace vanishes.
Contributions with $j\neq k$ in this regime can then only start at
$\mathcal{O}\left(\alpha_{s}^{2}(\overline{Q}_{3})\right)$ (and only through
$D=4$ field strength contributions, because of charge conjugation) and can
thus be neglected. We may thus write Eq. (3) as

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$
$\displaystyle=\sum_{j}\frac{ie_{q_{j}}^{2}}{e^{2}}\int\frac{d^{4}q_{4}}{(2\pi)^{4}}\int
d^{4}x_{1}\int d^{4}x_{2}\,e^{-i(q_{1}x_{1}+q_{2}x_{2})}$
$\displaystyle\qquad\qquad\times\langle
0|T(J_{j}^{\mu_{1}}(x_{1})J_{j}^{\mu_{2}}(x_{2}))|\gamma^{\mu_{3}}(q_{3})\gamma^{\mu_{4}}(q_{4})\rangle$
$\displaystyle\equiv\sum_{j}\frac{ie_{q_{j}}^{2}}{e^{2}}\Bigg{\langle}e^{-i(q_{1}x_{1}+q_{2}x_{2})}\,\bar{q}(x_{1})\gamma^{\mu_{1}}q(x_{1})\;\bar{q}(x_{2})\gamma^{\mu_{2}}q(x_{2})\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}_{q_{4},\,x_{1},\,x_{2}}\,.$
(16)

In this notation the subscripts on the angular bracket on the right refers to
integration variables, $j$ in the superscript to the specific quark-flavor in
the summation, and $\mu_{3},\mu_{4}$ indicate the matrix element with photons
in the final state. Defining the massless quark propagator $S(x)$ and its
Fourier transform $S(k)$ through

$S(x)=\int\frac{d^{4}k}{(2\pi)^{4}}\,e^{-ik\cdot
x}\,S(k)=\int\frac{d^{4}k}{(2\pi)^{4}}\,e^{-ik\cdot
x}\,\frac{\not{k}}{k^{2}}\,,$ (17)

we may contract two of the quark lines in Eq. (3) to yield

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}=$
$\displaystyle=\sum_{j}\frac{ie_{q_{j}}^{2}}{e^{2}}\Bigg{\langle}ie^{-i(q_{1}x_{1}+q_{2}x_{2})}\Big{[}\bar{q}(x_{1})\gamma^{\mu_{1}}S(x_{1}-x_{2})\gamma^{\mu_{2}}q(x_{2})+\bar{q}(x_{2})\gamma^{\mu_{2}}S(x_{2}-x_{1})\gamma^{\mu_{1}}q(x_{1})\Big{]}\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}_{q_{4},\,x_{1},\,x_{2}}$
$\displaystyle=\sum_{j}\frac{ie_{q_{j}}^{2}}{e^{2}}\Bigg{\langle}ie^{-i[x_{1}(k+q_{1})+x_{2}(-k+q_{2})]}\Big{[}\bar{q}(x_{1})\gamma^{\mu_{1}}S(k)\gamma^{\mu_{2}}q(x_{2})-\bar{q}(x_{2})\gamma^{\mu_{2}}S(k)\gamma^{\mu_{1}}q(x_{1})\Big{]}\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}_{q_{4},\,k,\,x_{1},\,x_{2}}\,.$
(18)

The two terms in the angular brackets take the form
$\bar{q}(x_{1})\Gamma\,q(x_{2})$ where $\Gamma$ is a Dirac matrix. We next use
translational invariance to transform to the symmetric point $(x_{1}-x_{2})/2$
according to

$\displaystyle\langle
0|\,\bar{q}(x_{1})\Gamma\,q(x_{2})\,|\gamma^{\mu_{3}}(q_{3})\gamma^{\mu_{4}}(q_{4})\rangle=$
$\displaystyle\,\langle 0|\,e^{i\hat{P}\cdot
x_{1}}\bar{q}(0)\Gamma\,e^{i\hat{P}\cdot(x_{2}-x_{1})}q(0)e^{-i\hat{P}\cdot
x_{2}}\,|\gamma(q_{3})^{\mu_{3}}\gamma^{\mu_{4}}(q_{4})\rangle$
$\displaystyle=$
$\displaystyle\,e^{-i(q_{3}+q_{4})\frac{x_{1}+x_{2}}{2}}\langle
0|\,\bar{q}\left(\frac{x_{1}-x_{2}}{2}\right)\Gamma\,q\left(\frac{x_{2}-x_{1}}{2}\right)\,|\gamma(q_{3})\gamma(q_{4})\rangle\,,$
(19)

where $\hat{P}$ is a momentum operator. The benefit of this translation is
that it allows to avoid total derivative operators in the OPE. At this stage
we can perform the OPE for $(x_{1}-x_{2})/2$ small. Through dimension $D=4$
one has

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$
$\displaystyle=-\sum_{j}\frac{e_{q_{j}}^{2}}{e^{2}}\Bigg{\langle}\bar{q}(0)\left[\Gamma^{\mu_{1}\mu_{2}}(-\hat{q})-\Gamma^{\mu_{2}\mu_{1}}(-\hat{q})\right]q(0)\Bigg{\rangle}^{j,\mu_{3},\mu_{4}}$
$\displaystyle-\sum_{j}\frac{e_{q_{j}}^{2}}{e^{2}}\Bigg{\langle}e^{-ix_{1}\left[k+q_{1}-\frac{q_{3}+q_{4}}{2}\right]}\,e^{-ix_{2}\left[-k+q_{2}-\frac{q_{3}+q_{4}}{2}\right]}$
$\displaystyle\qquad\qquad\times\frac{(x_{2}-x_{1})_{\alpha}}{2}\,\bar{q}(0)\left[\overrightarrow{\partial}^{\alpha}-\overleftarrow{\partial}^{\alpha}\right]\left[\Gamma^{\mu_{1}\mu_{2}}(k)+\Gamma^{\mu_{2}\mu_{1}}(k)\right]q(0)\Bigg{\rangle}^{j,\mu_{3},\mu_{4}}_{q_{4},x_{1},x_{2},k}\,$
$\displaystyle=-\sum_{j}\frac{e_{q_{j}}^{2}}{e^{2}}\,\Bigg{\langle}\bar{q}(0)\left[\Gamma^{\mu_{1}\mu_{2}}(-\hat{q})-\Gamma^{\mu_{2}\mu_{1}}(-\hat{q})\right]q(0)\Bigg{\rangle}^{j,\mu_{3},\mu_{4}}$
$\displaystyle+\sum_{j}\frac{e_{q_{j}}^{2}}{e^{2}}\,\lim_{p_{A}\to
0}\left(\frac{-i\,\partial^{p_{A}}_{\alpha}}{2}\right)\Bigg{\langle}\bar{q}(0)\left[\overrightarrow{\partial}^{\alpha}-\overleftarrow{\partial}^{\alpha}\right]\big{[}\Gamma^{\mu_{1}\mu_{2}}(\hat{q}+p_{A})$
$\displaystyle\hskip
150.0pt+\Gamma^{\mu_{2}\mu_{1}}(\hat{q}+p_{A})\big{]}q(0)\Bigg{\rangle}^{j,\mu_{3},\mu_{4}}\,,$
(20)

where $\Gamma^{\mu_{1}\mu_{2}}(k)=\gamma^{\mu_{1}}S(k)\gamma^{\mu_{2}}$ and we
introduced a spurious momentum $p_{A}$ to remove the $x_{i,\,\alpha}$ on the
second line through $x_{i,\,\alpha}=-i\lim_{p_{A}\to
0}\partial_{p_{A},\,\alpha}e^{ix_{i}p_{A}}$. For $D=4$ contributions we will
work in the chiral limit. Using the relation

$\Gamma^{\mu_{1}\mu_{2}}(k)+\Gamma^{\mu_{2}\mu_{1}}(k)=\frac{2\,k_{\delta}\gamma_{\beta}}{k^{2}}(g^{\mu_{1}\delta}g^{\mu_{2}\beta}+g^{\mu_{2}\delta}g^{\mu_{1}\beta}-g^{\mu_{1}\mu_{2}}g^{\delta\beta})\,.$
(21)

one finds

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$
$\displaystyle\approx-\sum_{j}\frac{e_{q_{j}}^{2}}{e^{2}}\,\Bigg{\langle}\bar{q}(0)\left[\Gamma^{\mu_{1}\mu_{2}}(-\hat{q})-\Gamma^{\mu_{2}\mu_{1}}(-\hat{q})\right]q(0)\Bigg{\rangle}^{j,\mu_{3},\mu_{4}}$
$\displaystyle-\sum_{j}\frac{ie_{q_{j}}^{2}}{e^{2}\hat{q}^{2}}\,\left(g^{\mu_{1}}_{\delta}g^{\mu_{2}}_{\beta}+g^{\mu_{2}}_{\delta}g^{\mu_{1}}_{\beta}-g^{\mu_{1}\mu_{2}}g_{\delta\beta}\right)\left(g_{\alpha}^{\delta}-2\,\frac{\hat{q}^{\delta}\hat{q}_{\alpha}}{\hat{q}^{2}}\right)$
$\displaystyle\qquad\qquad\times\Bigg{\langle}\bar{q}(0)\left[\overrightarrow{D}^{\alpha}-\overleftarrow{D}^{\alpha}\right]\gamma^{\beta}q(0)\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}\,.$
(22)

In the last line we have promoted the partial derivative to a covariant one,
which at the studied order is equivalent. Notice that the $D=4$ contribution
on the second and third lines already is a power correction of the $D=3$ one,
which at the same time is in a sub-leading kinematic regime for
$a_{\mu}^{\textrm{HLbL}}$. Imposing gauge invariance in both sides of the
equality, one can factor out $q_{4}$ to yield

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$
$\displaystyle=-q_{4\nu_{4}}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4\mu_{4}}}\,,$ $\displaystyle\langle Q\rangle^{j,\mu_{3}\mu_{4}}$
$\displaystyle=-q_{4\nu_{4}}\partial^{\mu_{4}}_{q_{4}}\langle
Q\rangle^{j,\mu_{3}\nu_{4}}\,,$ (23)

where $Q$ is one of the operators in the matrix elements in Eq. (3). The
tensor remaining after the factorisation has a well defined $q_{4}\rightarrow
0$ limit, and is antisymmetric under $\mu_{4}\leftrightarrow\nu_{4}$ exchange.
Since we have factored out the heavy scale $\hat{q}$, the remaining matrix
elements in the soft limit,

$\lim_{q_{4}\rightarrow 0}\partial^{\nu_{4}}_{q_{4}}\langle
Q\rangle^{j,\mu_{3}\mu_{4}}\,,$ (24)

can only depend on $q_{3}$. An alternative way of looking at this result is to
take the diagrams of Fig. 2 and expand in $p_{1},p_{2}$ small with
$q_{1,2}=\pm\hat{q}-(1/2)\left(p_{1}+p_{2}\right)$.

Figure 2: The diagrams needed for the lowest order in $\alpha_{s}$ OPE. The
diagram with $q_{1}$ and $q_{2}$ interchanged must be added.

A second kind of $D=4$ contribution arises from field strength tensors, which
are hidden in the implicit Dyson series of Eq. (3). Denoting the interaction
part of the action by $S_{\textrm{i}nt}$, one has a quadratic term

$e^{iS_{\mathrm{int}}}=\ldots-\frac{1}{2}\int d^{4}x_{3}\,\int
d^{4}x_{4}\,\mathcal{L}_{\mathrm{int}}(x_{3})\cdot\mathcal{L}_{\mathrm{int}}(x_{4})+\ldots\,,$
(25)

where the ellipses denote non-quadratic terms in the interaction Lagrangian

$\displaystyle\mathcal{L}_{\textrm{int}}(z_{i})$
$\displaystyle=\sum_{l=1}^{3}\bar{q}^{l}(z_{i})X_{l,\mu}\gamma^{\mu_{i}}q^{l}(z_{i})\,.$
(26)

Here, we include both gauge fields in $X_{\mu}$ as

$\displaystyle X_{l,\mu}(z)\equiv
g_{s}B^{a}_{\mu}(z)\frac{\lambda_{a}^{C}}{2}+e\,e_{q_{l}}\,A_{\mu}(z)\,,$ (27)

with $\lambda_{a}^{C}$ a Gell-Mann color matrix and $B_{\mu}^{a}(z)$ and
$A_{\mu}(z)$ are, respectively, the gluon and photon fields. The field
strength contributions are obtained by Taylor expanding the fields and keeping
the first gauge invariant non-zero order,

$\displaystyle A_{\mu}(x)$
$\displaystyle=\frac{1}{2}\,x^{\nu}F_{\nu\mu}=\frac{i}{2}\,F_{\nu\mu}\lim_{q_{4}\rightarrow
0}\partial_{q_{4}}^{\nu}e^{-iq_{4}x}\,,$ (28) $\displaystyle B_{\mu}^{a}(x)$
$\displaystyle=\frac{1}{2}\,x^{\nu}G_{\nu\mu}^{a}=\frac{i}{2}\,G_{\nu\mu}^{a}\lim_{q_{4}\rightarrow
0}\partial_{q_{4}}^{\nu}e^{-iq_{4}x}\,.$ (29)

Including these interaction terms in the Dyson series and closing the quark
lines for the leading (one-loop) order contribution, the color trace is
factored out and the associated Wilson coefficient contributions to both field
strength tensors are trivially related to the quark loop (i.e. the one-loop
HLbL tensor with a quark running around in the loop) according to

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$
$\displaystyle=\sum_{j}\frac{e_{q_{j}}^{2}}{e^{2}}\Bigg{\langle}e^{2}e_{q,j}^{2}F_{\nu_{3}^{\prime}\mu_{3}^{\prime}}F_{\nu_{4}^{\prime}\mu_{4}^{\prime}}+\frac{1}{2N_{c}}g_{s}^{2}G^{a}_{\nu_{3}^{\prime}\mu_{3}^{\prime}}G^{a}_{\nu_{4}^{\prime}\mu_{4}^{\prime}}\Bigg{\rangle}^{j,\mu_{3}\mu_{4}}$
$\displaystyle\qquad\quad\times\left(\frac{-1}{8}\lim_{q_{3},q_{4}\rightarrow
0}\,\partial_{q_{3}}^{\nu_{3}^{\prime}}\,\partial_{q_{4}}^{\nu_{4}^{\prime}}\,\Pi_{\mathrm{ql,j}}^{\mu_{1}\mu_{2}\mu_{3}^{\prime}\mu_{4}^{\prime}}\right)\,.$
(30)

Here the quark loop is given by

$\Pi_{\mathrm{ql},j}^{\mu_{1}\mu_{2}\mu_{3}^{\prime}\mu_{4}^{\prime}}=\left(-i\int\frac{d^{4}q_{4}}{(2\pi)^{4}}\left(\prod_{i=1}^{4}\int
d^{4}x_{i}\,e^{-iq_{i}x_{i}}\right)\langle
0|T\left(\prod_{j=1}^{4}J_{j}^{\mu_{j}}(x_{j})\right)|0\rangle\right)_{\mathrm{quark\,loop}}\,,$
(31)

with all the quark fields in the currents contracted. This topology is
diagrammatically depicted in Fig. 3. One can then factorise $q_{4}$ from the
matrix element (cf. Eq. (3)) and take the soft limit to yield the derivative
tensor

$\displaystyle\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial q_{4}^{\mu_{4}}}$
$\displaystyle=\sum_{j}\frac{e_{q_{j}}^{2}}{8e^{2}}\,\lim_{q_{4}\rightarrow
0}\left[\partial^{\nu_{4}}_{q_{4}}\,\Bigg{\langle}e^{2}e_{q,j}^{2}F_{\nu_{3}^{\prime}\mu_{3}^{\prime}}F_{\nu_{4}^{\prime}\mu_{4}^{\prime}}+\frac{1}{2N_{c}}g_{s}^{2}G^{a}_{\nu_{3}^{\prime}\mu_{3}^{\prime}}G^{a}_{\nu_{4}^{\prime}\mu_{4}^{\prime}}\Bigg{\rangle}^{j,\mu_{3}\mu_{4}}\right]$
$\displaystyle\qquad\quad\times\lim_{q_{3},q_{4}\rightarrow
0}\,\partial_{q_{3}}^{\nu_{3}^{\prime}}\,\partial_{q_{4}}^{\nu_{4}^{\prime}}\,\Pi_{\mathrm{ql,j}}^{\mu_{1}\mu_{2}\mu_{3}^{\prime}\mu_{4}^{\prime}}\,.$
(32)

An alternative way of looking at this is to take the diagrams of Fig. 3 and
expand in $q_{3},q_{4}$ small with
$q_{1,2}=\pm\hat{q}-(1/2)\left(q_{3}+q_{4}\right)$.

Figure 3: The diagrams needed for the lowest order in $\alpha_{s}$ OPE for the
operators with field strength tensors. The diagrams are with the $q_{3},q_{4}$
legs photon or gluon and all permutations of the external legs must be taken
into account.

Also factorising $q_{4}$ for the quark contributions in Eq. 3 and adding the
result to the field strength piece in Eq. 3, we then find at leading order
in333Let us note that while the gluon operator contribution is technically
$\alpha_{s}$-suppressed, and if $Q_{3}^{2}\gg\Lambda_{\mathrm{QCD}}^{2}$ the
associated contribution is going to enter suppressed by
$\alpha_{s}^{2}(Q_{3})$, there is no hard scale flowing through the associated
gluon lines, and then the associated contributions do not come with any
$\mathcal{O}(\alpha_{s}(\overline{Q}_{3}))$ suppression.
$\alpha_{s}(\overline{Q}_{3})$ and at NLO in the $Q_{3}/\overline{Q}_{3}$
expansion the result

$\displaystyle\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4,\,\mu_{4}}}\approx\sum_{j}\frac{e_{q_{j}}^{2}}{e^{2}}\,\lim_{q_{4}\rightarrow
0}\partial_{\nu_{4}}^{q_{4}}\Bigg{\langle}\bar{q}(0)\left[\Gamma^{\mu_{1}\mu_{2}}(-\hat{q})-\Gamma^{\mu_{2}\mu_{1}}(-\hat{q})\right]q(0)\Bigg{\rangle}^{j,\mu_{3},\mu_{4}}$
$\displaystyle+\sum_{j}\frac{ie_{q_{j}}^{2}}{e^{2}\hat{q}^{2}}\left(g_{\mu_{1}\delta}g_{\mu_{2}\beta}+g_{\mu_{2}\delta}g_{\mu_{1}\beta}-g_{\mu_{1}\mu_{2}}g_{\delta\beta}\right)\left(g^{\alpha\delta}-2\frac{\hat{q}^{\delta}\hat{q}^{\alpha}}{\hat{q}^{2}}\right)$
$\displaystyle\qquad\qquad\qquad\times\lim_{q_{4}\rightarrow
0}\partial_{\nu_{4}}^{q_{4}}\Bigg{\langle}\bar{q}(0)\left[\overrightarrow{D}^{\alpha}-\overleftarrow{D}^{\alpha}\right]\gamma^{\beta}q(0)\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}$
$\displaystyle+\sum_{j}\frac{e_{q_{j}}^{2}}{8e^{2}}\lim_{q_{4}\rightarrow
0}\left[\partial^{\nu_{4}}_{q_{4}}\Bigg{\langle}e^{2}e_{q,j}^{2}F_{\nu_{3}^{\prime}\mu_{3}^{\prime}}F_{\nu_{4}^{\prime}\mu_{4}^{\prime}}+\frac{1}{2N_{c}}g_{s}^{2}G^{a}_{\nu_{3}^{\prime}\mu_{3}^{\prime}}G^{a}_{\nu_{4}^{\prime}\mu_{4}^{\prime}}\Bigg{\rangle}^{j,\mu_{3}\mu_{4}}\right]$
$\displaystyle\qquad\qquad\qquad\times\lim_{q_{3},q_{4}\rightarrow
0}\partial_{q_{3}}^{\nu_{3}^{\prime}}\,\partial_{q_{4}}^{\nu_{4}^{\prime}}\,\Pi_{\mathrm{ql,j}}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}\,.$
(33)

At this stage the dependence on the two independent momenta, $\hat{q}$ and
$q_{3}$, appears completely factorized. The dependence on $\hat{q}$ is in the
Wilson coefficients, the dependence on $q_{3}$ in the matrix elements of the
operators. The matrix elements derived below are pictorially shown as a
Feynman diagram in Fig. 4.

Figure 4: The diagrams needed for the lowest order in $\alpha_{s}$ matrix
elements of the quark operators. The diagrams are with the $q_{3},q_{4}$ legs
interchanged must also be taken into account. The crossed circle indicates the
insertion of one of the quark operators.

##  4 Perturbative regime

If $Q_{3}^{2}\gg\Lambda_{\mathrm{QCD}}^{2}$, the matrix elements of the
operators in Eq. (33) can be computed perturbatively. In this section we
describe the process in some detail, also to illustrate some of the general
features of the OPE that can be trivially generalized to the nonperturbative
case, such as the need of renormalization and the exact power counting at each
step of the calculation, from the Lorentz decomposed HLbL tensor to the final
$\hat{\Pi}_{i}$.

###  4.1 Dimension $D=3$

At dimension $D=3$ we have

$\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4,\,\mu_{4}}}=\sum_{j}\frac{e_{q_{j}}^{2}}{e^{2}}\lim_{q_{4}\rightarrow
0}\partial^{\nu_{4}}_{q_{4}}\,\Bigg{\langle}\bar{q}(0)\left[\Gamma^{\mu_{1}\mu_{2}}(-\hat{q})-\Gamma^{\mu_{2}\mu_{1}}(-\hat{q})\right]q(0)\Bigg{\rangle}^{j,\mu_{3},\mu_{4}}\,.$
(34)

Expanding two electromagnetic vertices from the action one has

$\displaystyle\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4,\,\mu_{4}}}=$
$\displaystyle=\sum_{j}N_{c}\,e_{q,j}^{4}\lim_{q_{4}\rightarrow
0}i\,\partial_{q_{4}}^{\nu_{4}}\int\frac{d^{d}p}{(2\pi)^{d}}\,\mathrm{Tr}\left[\Gamma^{\mu_{1}\mu_{2}}(\hat{q})S(p+q_{3}+q_{4})\gamma^{\mu_{3}}S(p+q_{4})\gamma^{\mu_{4}}S(p)+(3\leftrightarrow
4)\right]\,,$ (35)

We may here take the derivative before integrating and use that for an
arbitrary momentum $k$ the propagator satisfies

$\partial_{q_{4}}^{\nu_{4}}S(k+q_{4})=-S(k+q_{4})\gamma^{\nu_{4}}S(k+q_{4})\,.$
(36)

In the $q_{4}\rightarrow 0$ limit, the remaining one-loop integral depending
on only the scale $q_{3}$ is then straightforward to compute. Alternatively
one can evaluate

$\lim_{q_{4}\rightarrow 0}\frac{\partial}{\partial q_{4,\,\mu_{4}}}\langle
0|\bar{q}(0)\left(\gamma^{\mu_{1}}\gamma^{\alpha}\gamma^{\mu_{2}}-\gamma^{\mu_{2}}\gamma^{\alpha}\gamma^{\mu_{1}}\right)q(0)|\gamma^{\mu_{3}}(q_{3})\gamma^{\nu_{4}}(q_{4})\rangle\,,$
(37)

via the Feynman diagrams of Fig. 4 and read off the Wilson coefficient needed
from Eq. (4.1).

At this level all the terms go as $\mathcal{O}(1/\hat{Q})$. A complication
arises when extracting the corresponding $\hat{\Pi}_{j}$ needed for
$a_{\mu}^{\textrm{HLbL}}$ in the integration in Eq. (10). The procedure
followed in Ref. [3] consisted in first extracting an intermediate set of $19$
independent scalar structures $\tilde{\Pi}_{j}$ in the tensor using
corresponding projectors $P^{\tilde{\Pi}_{j}}$. These $\tilde{\Pi}_{j}$ can
then be mapped onto the sought $\hat{\Pi}_{j}$ as in Ref. [3]. Examples of the
projectors to the $\tilde{\Pi}_{j}$ are
$P^{\tilde{\Pi}_{1}}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\nu_{4}}=g_{\mu_{1}\mu_{2}}\,g_{\mu_{3}\nu_{4}}\;q_{1,\mu_{4}}\,$
and
$P^{\tilde{\Pi}_{4}}_{\mu_{1}\mu_{2}\mu_{3}\mu_{4}\nu_{4}}=g_{\mu_{1}\mu_{2}}\,g_{\mu_{3}\nu_{4}}\;q_{2,\mu_{4}}$
, with the rest given in Ref. [3]. From these projectors, it is clear up to
which order one can unambiguously predict each $\tilde{\Pi}_{j}$: since
$q_{1}=\hat{q}-\frac{q_{3}}{2}$ and $q_{2}=-\hat{q}-\frac{q_{3}}{2}$, it is
easy to e.g. see that the dimension $D=3$ contribution is guaranteed to give
the correct answer for $\tilde{\Pi}_{1}$ and $\tilde{\Pi}_{4}$ up to and
including order $\left(\frac{1}{\overline{Q}_{3}}\right)^{0}$. Obviously, a
linear combination of $\tilde{\Pi}_{j}$ is unambiguously predicted at least up
to the most limiting of the addends. However further unambiguous terms can be
hidden through cancellation in linear combinations, which become lost if the
$\tilde{\Pi}_{j}$ have already been truncated. A simple example is
$\tilde{\Pi}_{1}+\tilde{\Pi}_{4}$, which is unambiguously predicted up to
order $\mathcal{O}\left(\frac{1}{\overline{Q}_{3}}\right)$.

In order to avoid the problem, we build a new set of $\tilde{\Pi}_{q_{3},\,j}$
with $q_{3}$ and $\hat{q}$ as variables, thus assuring we have a set with as
low possible powers of $\hat{Q}$ in the projectors. The advantage of this
basis is that the linear combination of $\tilde{\Pi}_{q_{3},\,j}$ is
unambiguously predicted at up to the most limiting of the addends, without
missing possible kinematic cancellations in the sum (no combination of
$\hat{q}$ and $q_{3}$ can lower the naive
$\mathcal{O}\left(\frac{1}{\overline{Q}_{3}}\right)$ order). The new
projectors to this basis are

$\displaystyle P^{\tilde{\Pi}_{1}}_{q_{3}}$
$\displaystyle=g_{\mu_{1}\mu_{2}}\,g_{\mu_{3}\nu_{4}}\,\hat{q}_{\mu_{4}}\,$
(38) $\displaystyle P^{\tilde{\Pi}_{2}}_{q_{3}}$
$\displaystyle=g_{\mu_{1}\mu_{2}}\,g_{\mu_{3}\nu_{4}}q_{3\mu_{4}}\,,$ (39)
$\displaystyle P^{\tilde{\Pi}_{3}}_{q_{3}}$
$\displaystyle=g_{\mu_{1}\mu_{3}}\,g_{\mu_{2}\nu_{4}}\,\hat{q}_{\mu_{4}}\,,$
(40) $\displaystyle P^{\tilde{\Pi}_{4}}_{q_{3}}$
$\displaystyle=g_{\mu_{1}\mu_{3}}\,g_{\mu_{2}\nu_{4}}q_{3\mu_{4}}\,,$ (41)
$\displaystyle P^{\tilde{\Pi}_{5}}_{q_{3}}$
$\displaystyle=g_{\mu_{2}\mu_{3}}\,g_{\mu_{1}\nu_{4}}\,\hat{q}_{\mu_{4}}\,,$
(42) $\displaystyle P^{\tilde{\Pi}_{6}}_{q_{3}}$
$\displaystyle=g_{\mu_{2}\mu_{3}}\,g_{\mu_{1}\nu_{4}}q_{3\mu_{4}}\,,$ (43)
$\displaystyle P^{\tilde{\Pi}_{7}}_{q_{3}}$
$\displaystyle=g_{\mu_{1}\nu_{4}}\,g_{\mu_{2}\mu_{4}}\,\hat{q}_{\mu_{3}}\,,$
(44) $\displaystyle P^{\tilde{\Pi}_{8}}_{q_{3}}$
$\displaystyle=g_{\mu_{2}\nu_{4}}\,g_{\mu_{3}\mu_{4}}q_{3\mu_{1}}\,,$ (45)
$\displaystyle P^{\tilde{\Pi}_{9}}_{q_{3}}$
$\displaystyle=g_{\mu_{3}\nu_{4}}\,g_{\mu_{1}\mu_{4}}q_{3\mu_{2}}\,,$ (46)
$\displaystyle P^{\tilde{\Pi}_{10}}_{q_{3}}$
$\displaystyle=g_{\mu_{1}\mu_{2}}\,\hat{q}_{\mu_{3}}\,\hat{q}_{\nu_{4}}q_{3\mu_{4}}\,,$
(47) $\displaystyle P^{\tilde{\Pi}_{11}}_{q_{3}}$
$\displaystyle=g_{\mu_{2}\mu_{3}}q_{3\mu_{1}}\,\hat{q}_{\nu_{4}}q_{3\mu_{4}}\,,$
(48) $\displaystyle P^{\tilde{\Pi}_{12}}_{q_{3}}$
$\displaystyle=g_{\mu_{3}\mu_{1}}q_{3\mu_{2}}q_{3\nu_{4}}\,\hat{q}_{\mu_{4}}\,,$
(49) $\displaystyle P^{\tilde{\Pi}_{13}}_{q_{3}}$
$\displaystyle=g_{\mu_{1}\nu_{4}}q_{3\mu_{2}}\,\hat{q}_{\mu_{3}}\,\hat{q}_{\mu_{4}}\,,$
(50) $\displaystyle P^{\tilde{\Pi}_{14}}_{q_{3}}$
$\displaystyle=g_{\mu_{1}\nu_{4}}q_{3\mu_{2}}\,\hat{q}_{\mu_{3}}q_{3\mu_{4}}\,,$
(51) $\displaystyle P^{\tilde{\Pi}_{15}}_{q_{3}}$
$\displaystyle=g_{\mu_{2}\nu_{4}}q_{3\mu_{1}}\,\hat{q}_{\mu_{3}}\,\hat{q}_{\mu_{4}}\,,$
(52) $\displaystyle P^{\tilde{\Pi}_{16}}_{q_{3}}$
$\displaystyle=g_{\mu_{2}\nu_{4}}q_{3\mu_{1}}\,\hat{q}_{\mu_{3}}q_{3\mu_{4}}\,,$
(53) $\displaystyle P^{\tilde{\Pi}_{17}}_{q_{3}}$
$\displaystyle=g_{\mu_{3}\nu_{4}}q_{3\mu_{1}}q_{3\mu_{2}}\,\hat{q}_{\mu_{4}}\,,$
(54) $\displaystyle P^{\tilde{\Pi}_{18}}_{q_{3}}$
$\displaystyle=g_{\mu_{3}\nu_{4}}q_{3\mu_{1}}q_{3\mu_{2}}q_{3\mu_{4}}\,,$ (55)
$\displaystyle P^{\tilde{\Pi}_{19}}_{q_{3}}$
$\displaystyle=q_{3\mu_{1}}q_{3\mu_{2}}\,\hat{q}_{\mu_{3}}\,\hat{q}_{\nu_{4}}q_{3\mu_{4}}\,,$
(56)

from which relations can be found between the old $\tilde{\Pi}_{j}$ and the
new $\tilde{\Pi}_{q_{3},\,j}$, as needed to reach the $\hat{\Pi}_{j}$. We
repeat exactly the same procedure for the case in which $q_{1}$ or $q_{2}$ is
the small momentum. Notice how building the new $q_{1}$ or $q_{2}$ projectors
can be done by cyclic permutation, i.e. by changing $q_{3}\rightarrow q_{1}$
and
$\hat{q}=\left(q_{1}-q_{2}\right)/2\rightarrow\hat{q}_{23}=\left(q_{2}-q_{3}\right)/2$
for the $q_{1}$ small case, and $q_{3}\rightarrow q_{2}$ and
$\hat{q}=\left(q_{1}-q_{2}\right)/2\rightarrow\hat{q}_{31}=\left(q_{3}-q_{1}\right)/2$
for the $q_{2}$ small case.

At this stage it is straightforward to extract the linear combinations of
$\tilde{\Pi}_{q_{3},\,j}$ that can enter order by order in the
$\hat{\Pi}_{j}$. We find that truncating at order $D$ in the operators one can
unambiguously predict,

$\displaystyle\hat{\Pi}_{1}\rightarrow\frac{1}{\overline{Q}_{3}^{D-1}}\,,\quad\hat{\Pi}_{4}\rightarrow\frac{1}{\overline{Q}_{3}^{D}}\,,\quad\hat{\Pi}_{7}\rightarrow\frac{1}{\overline{Q}_{3}^{D+1}}\,,\quad\hat{\Pi}_{17}\rightarrow\frac{1}{\overline{Q}_{3}^{D}}\,,$
$\displaystyle\hat{\Pi}_{39}\rightarrow\frac{1}{\overline{Q}_{3}^{D}}\,,\quad\hat{\Pi}_{54}\rightarrow\frac{1}{\overline{Q}_{3}^{D}}\,.$
(57)

At $D=3$ one then finds for $q_{3}$ small

$\displaystyle\underline{q_{3}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=-\frac{4N_{c}\sum_{j}e_{q,j}^{4}}{\pi^{2}Q_{3}^{2}\,\overline{Q}_{3}^{2}}+\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{3}}\right)\,,\,\hat{\Pi}_{4}=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{4}}\right)\,,\,\hat{\Pi}_{7}=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{5}}\right)\,,\,$
(58) $\displaystyle\hat{\Pi}_{17}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{4}}\right)\,,\,\hat{\Pi}_{39}=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{4}}\right)\,,\,\hat{\Pi}_{54}=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{4}}\right)\,.$
(59)

For $q_{1}$ small the result is

$\displaystyle\underline{q_{1}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{4}}\right)\,,\,\hat{\Pi}_{4}=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{2}}\right)\,,\,\hat{\Pi}_{7}=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{3}}\right)\,,\,$
(60) $\displaystyle\hat{\Pi}_{17}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{4}}\right)\,,\,\hat{\Pi}_{39}=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{4}}\right)\,,\,\hat{\Pi}_{54}=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{4}}\right)\,.$
(61)

Finally, for $q_{2}$ small we find

$\displaystyle\underline{q_{2}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{4}}\right)\,,\,\hat{\Pi}_{4}=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{2}}\right)\,,\,\hat{\Pi}_{7}=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{4}}\right)\,,\,$
(62) $\displaystyle\hat{\Pi}_{17}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{4}}\right)\,,\,\hat{\Pi}_{39}=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{4}}\right)\,,\,\hat{\Pi}_{54}=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{4}}\right)\,.$
(63)

###  4.2 Dimension $D=4$

The first operator appearing at dimension $D=4$ in Eq. (33) involves a quark
current. Expanding down two electromagnetic vertices from the Dyson series in
the matrix element and contracting one finds, in the perturbative limit
$Q_{3}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$,

$\displaystyle\lim_{q_{4}\rightarrow
0}\partial^{\nu_{4}}_{q_{4}}\,\Bigg{\langle}\bar{q}(0)\left[\overrightarrow{D}^{\alpha}-\overleftarrow{D}^{\alpha}\right]\gamma^{\beta}q(0)\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}$
$\displaystyle=-2N_{c}\,e^{2}e_{q,j}^{2}\,\lim_{q_{4}\rightarrow
0}\,\partial^{\nu_{4}}_{q_{4}}\int\frac{d^{4}k}{(2\pi)^{4}}\,k^{\alpha}\,\mathrm{Tr}[\gamma^{\beta}S(k+q_{3}+q_{4})\gamma^{\mu_{3}}S(k+q_{4})\gamma^{\mu_{4}}S(k)+(3\leftrightarrow
4)]\,.$ (64)

The result for the matrix element is divergent at $d=4$. We regularize the
divergence within dimensional regularization, with $d=4-2\epsilon$. Inserting
Eq. (64) into Eq. (33), one finds the derivative of the HLbL tensor, where the
divergence in the matrix element in Eq. (64) cancels the divergence from the
Wilson coefficient of the field strength contribution at dimension $D=4$ in
Eq. 33, thus yielding a finite result.

An alternative approach is to separately project the quark operator matrix
element and the photon matrix element on the $\hat{\Pi}_{i}$. As we will
explain below, one also has to re-expand the $D=3$ contribution in Eq. (33).
In the following we study the three OPE contributions separately.

####  4.2.1 Dimension $D=4$ quark current operator

The $\hat{\Pi}_{i}$ are finite in the $q_{3}$-small regime. Leaving out the
overall factor $N_{c}\sum_{j}e_{q_{j}}^{4}$ for the $\hat{\Pi}_{i}$ here and
in the following, one finds

$\displaystyle\underline{q_{3}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{4}}\right)\,,\qquad\qquad\qquad\hat{\Pi}_{4\,}=\frac{16(\delta_{12}^{4}-\delta_{12}^{2}Q_{3}^{2}+Q_{3}^{4})}{3\pi^{2}(Q_{3}^{2}-\delta_{12}^{2})^{2}\overline{Q}_{3}^{4}}\,,\,$
$\displaystyle\hat{\Pi}_{7}$
$\displaystyle=\frac{32\delta_{12}Q_{3}^{2}}{3\pi^{2}(Q_{3}^{2}-\delta_{12}^{2})^{2}\overline{Q}_{3}^{5}}\,,\,\qquad\hat{\Pi}_{17}=\frac{4(Q_{3}^{2}-4\delta_{12}^{2})}{3\pi^{2}Q_{3}^{2}(Q_{3}^{2}-\delta_{12}^{2})\overline{Q}_{3}^{4}}\,,\,$
$\displaystyle\hat{\Pi}_{39}$
$\displaystyle=\frac{4(Q_{3}^{4}+4\delta_{12}^{4}-9\delta_{12}^{2}Q_{3}^{2})}{3\pi^{2}Q_{3}^{2}(Q_{3}^{2}-\delta_{12}^{2})^{2}\overline{Q}_{3}^{4}}\,,\,\hat{\Pi}_{54}=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{5}}\right)\,.$
(65)

It is important to note that there are kinematic singularities in the limit
$\delta_{12}\rightarrow Q_{3}$. Since the $\hat{\Pi}_{i}$ by construction
should be free of such singularities, they must cancel against the remaining
contribution from the OPE in Eq. (33), thus connecting the different $D=3$ and
$D=4$ Green functions. In this section we will show how this explicitly occurs
in the perturbative regime. The fact that it must also occur in the
nonperturbative regime is going to lead to simple nontrivial relations among
the different form factors appearing in Lorentz decompositions of the
nonperturbative matrix elements, as explained in Sec. 5.

An additional feature for the $\hat{\Pi}_{i}$ appears when $q_{1}$ or $q_{2}$
is small. With $1/\hat{\epsilon}=1/\epsilon-\gamma_{\textrm{E}}+\log 4\pi$ and
neglecting terms of order $\epsilon$ and higher, one finds

$\displaystyle\underline{q_{1}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=-\frac{16[3+\epsilon(-4-6\ln
Q_{1})]}{9\pi^{2}\hat{\epsilon}\,\overline{Q}_{1}^{4}}\,,\,\hat{\Pi}_{4}=-\frac{4\delta_{23}^{4}-9\delta_{23}^{2}Q_{1}^{2}+Q_{1}^{4}}{3\pi^{2}Q_{1}^{2}(Q_{1}^{2}-\delta_{23}^{2})^{2}\overline{Q}_{1}^{2}}\,,\,$
$\displaystyle\hat{\Pi}_{7}$
$\displaystyle=4\frac{3\delta_{23}^{3}-5\delta_{23}Q_{1}^{2}}{3\pi^{2}Q_{1}^{2}(Q_{1}^{2}-\delta_{23}^{2})^{2}\overline{Q}_{1}^{3}}\,,\,\quad\hat{\Pi}_{17}=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{5}}\right)\,,\,$
$\displaystyle\hat{\Pi}_{39}$
$\displaystyle=4\frac{4\delta_{23}^{4}-9\delta_{23}^{2}Q_{1}^{2}+Q_{1}^{4}}{3\pi^{2}Q_{1}^{2}(Q_{1}^{2}-\delta_{23}^{2})^{2}\overline{Q}_{1}^{4}}\,,\,\quad\hat{\Pi}_{54}=-4\frac{Q_{1}^{2}+2\delta_{23}^{2}}{3\pi^{2}Q_{1}^{2}(Q_{1}^{2}-\delta_{23}^{2})\overline{Q}_{1}^{4}}\,.$
(66) $\displaystyle\underline{q_{2}\textrm{ small:}}$
$\displaystyle\hat{\Pi}_{1}$ $\displaystyle=-\frac{16[3+\epsilon(-4-6\ln
Q_{2})]}{9\pi^{2}\hat{\epsilon}\,\overline{Q}_{2}^{4}}\,,\,\hat{\Pi}_{4}=-\frac{4\delta_{31}^{4}-9\delta_{31}^{2}Q_{2}^{2}+Q_{2}^{4}}{3\pi^{2}Q_{2}^{2}(Q_{2}^{2}-\delta_{31}^{2})^{2}\overline{Q}_{2}^{2}}\,,\,$
$\displaystyle\hat{\Pi}_{7}$
$\displaystyle=-4\frac{4\delta_{31}^{4}-9\delta_{31}^{2}Q_{2}^{2}+Q_{2}^{4}}{3\pi^{2}Q_{2}^{2}(Q_{2}^{2}-\delta_{31}^{2})^{2}\overline{Q}_{2}^{4}}\,,\,\hat{\Pi}_{17}=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{5}}\right)\,,\,$
$\displaystyle\hat{\Pi}_{39}$
$\displaystyle=4\frac{4\delta_{31}^{4}-9\delta_{31}^{2}Q_{2}^{2}+Q_{2}^{4}}{3\pi^{2}Q_{2}^{2}(Q_{2}^{2}-\delta_{31}^{2})^{2}\overline{Q}_{2}^{4}}\,,\,\quad\hat{\Pi}_{54}=4\frac{Q_{2}^{2}+2\delta_{31}^{2}}{3\pi^{2}Q_{2}^{2}(Q_{2}^{2}-\delta_{31}^{2})\overline{Q}_{2}^{4}}\,.$
(67)

As can be seen, $\hat{\Pi}_{1}$ contains a UV divergent logarithm in both
limits.

####  4.2.2 Dimension $D=4$ photon field strength operator

Taking into account that at the studied order no other matrix elements are
divergent, the only way of cancelling the logarithmic divergences in Eqs.
(4.2.1)–(4.2.1) is with a similar (in this case infrared) logarithm in a
Wilson coefficient. This comes from the photon field strength contribution,

$\displaystyle\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4,\,\mu_{4}}}=\frac{\sum_{j}e_{q_{j}}^{4}}{8}\,\lim_{q_{4}\rightarrow
0}\left[\partial^{\nu_{4}}_{q_{4}}\,\left\langle
F_{\nu_{3}^{\prime}\mu_{3}^{\prime}}F_{\nu_{4}^{\prime}\mu_{4}^{\prime}}\right\rangle^{\mu_{3}\mu_{4}}\right]\,\lim_{q_{3},q_{4}\rightarrow
0}\partial_{q_{3}}^{\nu_{3}^{\prime}}\,\partial_{q_{4}}^{\nu_{4}^{\prime}}\,\Pi_{\mathrm{ql}}^{\mu_{1}\mu_{2}\mu_{3}^{\prime}\mu_{4}^{\prime}}\,.$

Since photons do not see strong interactions at the studied order in $\alpha$,
the computation of the associated matrix element is straightforward. The
result is

$\displaystyle\lim_{q_{4}\rightarrow
0}\partial^{\nu_{4}}_{q_{4}}\,\left\langle
F_{\nu_{3}^{\prime}\mu_{3}^{\prime}}F_{\nu_{4}^{\prime}\mu_{4}^{\prime}}\right\rangle^{\mu_{3}\mu_{4}}=$
$\displaystyle-\Bigg{[}\left(q_{3,\,\mu_{3}^{\prime}}\,g_{\nu_{3}^{\prime}}^{\mu_{3}}-q_{3,\,\nu_{3}^{\prime}}\,g_{\mu_{3}^{\prime}}^{\mu_{3}}\right)\left(g_{\mu_{4}^{\prime}}^{\nu_{4}}\,g_{\nu_{4}^{\prime}}^{\mu_{4}}-g_{\nu_{4}^{\prime}}^{\nu_{4}}\,g_{\mu_{4}^{\prime}}^{\mu_{4}}\right)$
$\displaystyle\quad+\left(g_{\mu_{3}^{\prime}}^{\nu_{4}}\,g_{\nu_{3}^{\prime}}^{\mu_{4}}-g_{\nu_{3}^{\prime}}^{\nu_{4}}\,g_{\mu_{3}^{\prime}}^{\mu_{4}}\right)\left(q_{3,\,\mu_{4}^{\prime}}\,g_{\nu_{4}^{\prime}}^{\mu_{3}}-q_{3,\,\nu_{4}^{\prime}}\,g_{\mu_{4}^{\prime}}^{\mu_{3}}\right)\Bigg{]}\,.$
(68)

The difficulty here lies in the computation of

$\displaystyle\lim_{q_{3},q_{4}\rightarrow
0}\partial_{q_{3}}^{\nu_{3}^{\prime}}\,\partial_{q_{4}}^{\nu_{4}^{\prime}}\,\Pi_{\mathrm{ql}}^{\mu_{1}\mu_{2}\mu_{3}^{\prime}\mu_{4}^{\prime}}\,.$
(69)

The massless quark loop tensor,
$\Pi_{\mathrm{ql}}^{\mu_{1}\mu_{2}\mu_{3}^{\prime}\mu_{4}^{\prime}}$ is well-
known to be finite and also the derivative $\lim_{q_{4}\rightarrow
0}\partial_{q_{4}}^{\nu_{4}^{\prime}}\,\Pi_{\mathrm{ql}}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}^{\prime}}$.
However, if one tries to apply the remaining $\lim_{q_{3}\rightarrow
0}\partial_{q_{3}}^{\nu_{3}^{\prime}}$ operation to the known analytic
expression of $\lim_{q_{4}\rightarrow
0}\partial_{q_{4}}^{\nu_{4}^{\prime}}\Pi_{\mathrm{ql}}^{\mu_{1}\mu_{2}\mu_{3}^{\prime}\mu_{4}^{\prime}}$,
one finds divergent logarithms. Indeed, these IR divergences at the Wilson
coefficient level are the ones that must cancel the corresponding UV
divergences of the quark current matrix element discussed above. To recover
the correct finite terms, however, one needs to regularize both the Wilson
coefficient and the matrix element in the same way, which we do using
dimensional regularization. Instead of first performing the multi-scale loop
integral and then applying the derivatives and limits $\lim_{q_{i}\rightarrow
0}\partial_{q_{i}}^{\nu_{i}^{\prime}}$, we once again first apply the
operations and then straightforwardly calculate the resulting single-scale
massless one-loop integrals. Combining them together one finds the associated
derivative HLbL tensor. Projecting them into the $\hat{\Pi}_{i}$ one
respectively obtains in the three regions of small $q_{i}$

$\displaystyle\underline{q_{3}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{4}}\right)\,,\,\qquad\qquad\quad\hat{\Pi}_{4}=-\frac{16(2\delta_{12}^{4}-3\delta_{12}^{2}Q_{3}^{2}+2Q_{3}^{4})}{3\pi^{2}(Q_{3}^{2}-\delta_{12}^{2})^{2}\overline{Q}_{3}^{4}}\,,\,$
$\displaystyle\hat{\Pi}_{7}$
$\displaystyle=-\frac{32\delta_{12}Q_{3}^{2}}{3\pi^{2}(Q_{3}^{2}-\delta_{12}^{2})^{2}\overline{Q}_{3}^{5}}\,,\,\hat{\Pi}_{17}=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{5}}\right)\,,\,$
$\displaystyle\hat{\Pi}_{39}$
$\displaystyle=\frac{16\delta_{12}^{2}}{3\pi^{2}(Q_{3}^{2}-\delta_{12}^{2})^{2}\overline{Q}_{3}^{4}}\,,\,\quad\hat{\Pi}_{54}=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{5}}\right)\,,$
(70) $\displaystyle\underline{q_{1}\textrm{ small:}}$
$\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\frac{16[1+\epsilon(-3-2\ln\frac{\overline{Q}_{1}}{2})]}{3\pi^{2}\hat{\epsilon}\,\overline{Q}_{1}^{4}}\,,\,\quad\hat{\Pi}_{4}=-\frac{4\delta_{23}^{2}}{3\pi^{2}(Q_{1}^{2}-\delta_{23}^{2})^{2}\overline{Q}_{1}^{2}}\,,\,$
$\displaystyle\hat{\Pi}_{7}$
$\displaystyle=8\frac{\delta_{23}}{3\pi^{2}(Q_{1}^{2}-\delta_{23}^{2})^{2}\overline{Q}_{1}^{3}}\,,\,\qquad\hat{\Pi}_{17}=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{5}}\right)\,,\,$
$\displaystyle\hat{\Pi}_{39}$
$\displaystyle=\frac{16\,\delta_{23}^{2}}{3\pi^{2}(Q_{1}^{2}-\delta_{23}^{2})^{2}\overline{Q}_{1}^{4}}\,,\,\qquad\quad\hat{\Pi}_{54}=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{5}}\right)\,,$
(71) $\displaystyle\underline{q_{2}\textrm{ small:}}$
$\displaystyle\hat{\Pi}_{1}$
$\displaystyle=16\frac{1+\epsilon(-3-2\ln\frac{\overline{Q}_{2}}{2})}{3\pi^{2}\hat{\epsilon}\,\overline{Q}_{2}^{4}}\,,\,\quad\hat{\Pi}_{4}=-\frac{4\delta_{31}^{2}}{3\pi^{2}(Q_{2}^{2}-\delta_{31}^{2})^{2}\overline{Q}_{2}^{2}}\,,\,$
$\displaystyle\hat{\Pi}_{7}$
$\displaystyle=-\frac{16\,\delta_{31}^{2}}{3\pi^{2}(Q_{2}^{2}-\delta_{31}^{2})^{2}\overline{Q}_{2}^{2}}\,,\,\quad\hat{\Pi}_{17}=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{5}}\right)\,,\,$
$\displaystyle\hat{\Pi}_{39}$
$\displaystyle=\frac{16\,\delta_{31}^{2}}{3\pi^{2}(Q_{2}^{2}-\delta_{31}^{2})^{2}\overline{Q}_{2}^{4}}\,,\,\quad\quad\hat{\Pi}_{54}=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{5}}\right)\,.$
(72)

Comparing to the $\hat{\Pi}_{i}$ for the quark-current matrix element in Eqs.
(4.2.1)–(4.2.1), one notices that the divergence in $\epsilon$cancels when
adding the contributions fro the quark and field strength operators. The
infrared and ultraviolet logarithms nicely recombine into
$\log\frac{\overline{Q}_{i}}{2Q_{i}}$. However, the collinear divergences when
$\delta_{ij}\rightarrow Q_{k}$ remain.

####  4.2.3 Expanded dimension $D=3$ operator

To cancel the collinear divergences, one also has to add the previously
computed $D=3$ expanded to NLO in the energy expansion for the
$\hat{\Pi}_{i}$. The $\hat{\Pi}_{i}$ in the three momentum regimes, including
also the previous LO $\hat{\Pi}_{1}$ contribution from Eq. (4.1), are given by

$\displaystyle\underline{q_{3}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=-\frac{4}{\pi^{2}\,Q_{3}^{2}\,\overline{Q}_{3}^{2}}\,,\,\quad\hat{\Pi}_{4}=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{5}}\right)\,,\,\quad\hat{\Pi}_{7}=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{6}}\right)\,,\,$
$\displaystyle\hat{\Pi}_{17}$
$\displaystyle=\frac{4}{\pi^{2}(Q_{3}^{2}-\delta_{12}^{2})\overline{Q}_{3}^{4}}\,,\,\quad\hat{\Pi}_{39}=\frac{4}{\pi^{2}(Q_{3}^{2}-\delta_{12}^{2})\overline{Q}_{3}^{4}}\,,\,\quad\hat{\Pi}_{54}=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{5}}\right)\,.$
(73) $\displaystyle\underline{q_{1}\textrm{ small:}}$
$\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{5}}\right)\,,\,\quad\hat{\Pi}_{4}=-\frac{1}{\pi^{2}(Q_{1}^{2}-\delta_{23}^{2})\overline{Q}_{1}^{2}}\,,\,\quad\hat{\Pi}_{7}=\frac{4\delta_{23}}{\pi^{2}(Q_{1}^{2}-\delta_{23}^{2})Q_{1}^{2}\overline{Q}_{1}^{3}}\,,\,$
$\displaystyle\hat{\Pi}_{17}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{5}}\right)\,,\,\quad\hat{\Pi}_{39}=\frac{4}{\pi^{2}(Q_{1}^{2}-\delta_{23})\overline{Q}_{1}^{4}}\,,\,\quad\hat{\Pi}_{54}=\frac{4}{\pi^{2}(Q_{1}^{2}-\delta_{23}^{2})\overline{Q}_{1}^{4}}\,.$
(74) $\displaystyle\underline{q_{2}\textrm{ small:}}$
$\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{5}}\right)\,,\,\quad\hat{\Pi}_{17}=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{5}}\right)\,,\,\quad\hat{\Pi}_{4}=-\frac{1}{\pi^{2}(Q_{2}^{2}-\delta_{31}^{2})\overline{Q}_{2}^{2}}\,,\,$
$\displaystyle\hat{\Pi}_{7}$
$\displaystyle=-\hat{\Pi}_{39}=\hat{\Pi}_{54}=-\frac{4}{\pi^{2}(Q_{2}^{2}-\delta_{31}^{2})\overline{Q}_{2}^{4}}\,.$
(75)

####  4.2.4 Final $\hat{\Pi}_{i}$ at dimension $D=4$

At this stage we can add the three contributions together (Eqs. (4.2.1–4.2.2)
for the quark current operator, Eqs. (4.2.2–4.2.2) for the photon operator and
Eqs. (4.2.3–4.2.3)) in the three momentum limits to get a finite result at NLO
in the $\frac{1}{\overline{Q}_{3}}$ expansion. The results are

$\displaystyle\underline{q_{3}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=-\frac{4}{\pi^{2}\,Q_{3}^{2}\,\overline{Q}_{3}^{2}}\quad\hat{\Pi}_{4}=-\frac{16}{3\pi^{2}\overline{Q}_{3}^{4}}\,,\,\quad\hat{\Pi}_{7}=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{6}}\right)\,,\,$
$\displaystyle\hat{\Pi}_{17}$
$\displaystyle=\frac{16}{3\pi^{2}Q_{3}^{2}\overline{Q}_{3}^{4}}\,,\,\quad\hat{\Pi}_{39}=\frac{16}{3\pi^{2}Q_{3}^{2}\overline{Q}_{3}^{4}}\,,\,\quad\hat{\Pi}_{54}=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{5}}\right)\,,$
(76) $\displaystyle\underline{q_{1}\textrm{ small:}}$
$\displaystyle\hat{\Pi}_{1}$
$\displaystyle=-\frac{16\left(5+6\ln\frac{\overline{Q}_{1}}{2Q_{1}}\right)}{9\pi^{2}\overline{Q}_{1}^{4}}\,,\,\quad\hat{\Pi}_{4}=-\frac{4}{3\pi^{2}Q_{1}^{2}\overline{Q}_{1}^{2}}\,,\,$
$\displaystyle\hat{\Pi}_{7}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{4}}\right)\,,\,\quad\qquad\hat{\Pi}_{17}=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{5}}\right)\,,\,$
$\displaystyle\hat{\Pi}_{39}$
$\displaystyle=\frac{16}{3\pi^{2}Q_{1}^{2}\overline{Q}_{1}^{4}}\,,\,\quad\hat{\Pi}_{54}=\frac{8}{3\pi^{2}Q_{1}^{2}\overline{Q}_{1}^{4}}\,,$
(77) $\displaystyle\underline{q_{2}\textrm{ small:}}$
$\displaystyle\hat{\Pi}_{1}$
$\displaystyle=-\frac{16\left(5+6\ln\frac{\overline{Q}_{2}}{2Q_{2}}\right)}{9\pi^{2}\overline{Q}_{2}^{4}}\,,\,\quad\hat{\Pi}_{4}=-\frac{4}{3\pi^{2}Q_{2}^{2}\overline{Q}_{2}^{2}}\,,\,$
$\displaystyle\hat{\Pi}_{7}$
$\displaystyle=-\frac{16}{3\pi^{2}Q_{2}^{2}\overline{Q}_{2}^{4}}\,,\,\qquad\qquad\hat{\Pi}_{17}=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{5}}\right)\,,\,$
$\displaystyle\hat{\Pi}_{39}$
$\displaystyle=\frac{16}{3\pi^{2}Q_{2}^{2}\overline{Q}_{2}^{4}}\,,\,\qquad\qquad\quad\hat{\Pi}_{54}=-\frac{8}{3\pi^{2}Q_{2}^{2}\overline{Q}_{2}^{4}}\,.$
(78)

The OPE developed in Refs. [1, 2, 3] was valid in the purely short-distance
limit $Q_{i}^{2},Q_{j}^{2},Q_{k}^{2}\gg\Lambda_{\mathrm{QCD}}^{2}$. Here we
have considered the Melnikov-Vainshtein limit $Q_{i}^{2},Q_{j}^{2}\gg
Q_{k}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$ where also the momentum $Q_{k}^{2}$ is
in the perturbative regime. This is part of the region444Notice however how in
the very extreme (and phenomenologically irrelevant) case in which
$Q_{1,2}^{2}\gg Q_{3}^{2}\gg\Lambda_{\mathrm{QCD}}^{2}$, one may need to resum
$\ln\frac{\overline{Q}_{3}}{2Q_{3}}$ factors.
$Q_{i}^{2},Q_{j}^{2},Q_{k}^{2}\gg\Lambda_{\mathrm{QCD}}^{2}$ and was called
the corner region in Ref. [3]. The $\hat{\Pi}_{i}$ in Eqs. (4.2.4–4.2.4)
should therefore agree with the corner region $\hat{\Pi}_{i}$ in Ref. [3].
Indeed, all the $\hat{\Pi}_{i}$ exactly agree.

####  4.2.5 Renormalization of the OPE

Up to now we have worked with divergent matrix elements and Wilson
coefficients. One can instead work with both finite matrix elements and Wilson
coefficients by performing standard renormalization of the OPE, see e.g. Ref.
[43]. At first sight it may seem striking that already at the leading order
contribution in $\alpha_{s}$ and $\alpha$ one has to perform such a task for
the dimension $D=4$ contribution. This can be regarded as a consequence of an
intermediate mismatch in the order in $\alpha$ of the two contributing
operators, of the respective forms $\bar{q}D\Gamma q$ and $\alpha\,F\cdot F$
where $\Gamma$ is a Dirac matrix, which eventually enter at the same order in
$\alpha$, since the former needs to expand two extra EM vertices in order to
capture the two initial photons. In that intermediate step, the QED mixing
between $\bar{q}D\Gamma q$ and $\alpha\,F\cdot F$ needs to be considered, even
when the self-renormalization of each separate operator is not going to affect
our calculation at the studied order.

In order to assess the mixing, one can simply focus on the loop divergences of
the relevant matrix element. For the quark current operator we have

$\displaystyle\lim_{q_{4}\rightarrow
0}\partial^{\nu_{3}}_{q_{3}}\,\partial^{\nu_{4}}_{q_{4}}\,\Bigg{\langle}\bar{q}(0)\left[\overrightarrow{D}^{\alpha}-\overleftarrow{D}^{\alpha}\right]\gamma^{\beta}q(0)\Bigg{\rangle}^{\mu_{3},\,\mu_{4},j}={\textrm{finite}}-i\,\frac{8}{3}\,\frac{\alpha}{4\pi}\,N_{c}\,e_{q,j}^{2}\,\frac{Q_{3}^{-2\epsilon}}{\epsilon}$
$\displaystyle\times\Big{(}-g^{\mu_{3}\mu_{4}}g^{\alpha\beta}g^{\nu_{3}\nu_{4}}+g^{\mu_{3}\mu_{4}}g^{\alpha\nu_{3}}g^{\beta\nu_{4}}+g^{\mu_{3}\mu_{4}}g^{\alpha\nu_{4}}g^{\beta\nu_{3}}+g^{\mu_{3}\alpha}g^{\mu_{4}\beta}g^{\nu_{3}\nu_{4}}-g^{\mu_{3}\alpha}g^{\mu_{4}\nu_{3}}g^{\beta\nu_{4}}$
$\displaystyle+g^{\mu_{3}\beta}g^{\mu_{4}\alpha}g^{\nu_{3}\nu_{4}}-g^{\mu_{3}\beta}g^{\mu_{4}\nu_{3}}g^{\alpha\nu_{4}}-g^{\mu_{3}\nu_{4}}g^{\mu_{4}\alpha}g^{\beta\nu_{3}}-g^{\mu_{3}\nu_{4}}g^{\mu_{4}\beta}g^{\alpha\nu_{3}}+g^{\mu_{3}\nu_{4}}g^{\mu_{4}\nu_{3}}g^{\alpha\beta}\Big{)}$
$\displaystyle=\textrm{finite}-i\,\frac{2}{3}\,N_{c}\,e_{q,j}^{2}\,\frac{Q_{3}^{-2\epsilon}}{\hat{\epsilon}}\lim_{q_{4}\rightarrow
0}\,\partial^{\nu_{3}}_{q_{3}}\,\partial^{\nu_{4}}_{q_{4}}\,\Bigg{\langle}\frac{\alpha}{4\pi}\left(F^{\mu\nu}F_{\mu\nu}g^{\alpha\beta}+d\,F^{\alpha\gamma}F_{\gamma}^{\beta}\right)\Bigg{\rangle}^{\mu_{3}\mu_{4}}\,.$
(79)

In the final step we recognised the photon field strength matrix element. The
$\overline{\mathrm{MS}}$ prescription consists in removing the
$\frac{1}{\hat{\epsilon}}$ pole from the matrix element and incorporating a
scale $\mu$ to make the implicit logarithm in $Q_{3}^{-2\epsilon}$
dimensionless, i.e.

$\displaystyle\lim_{q_{4}\rightarrow
0}\partial^{\nu_{3}}_{q_{3}}\,\partial^{\nu_{4}}_{q_{4}}\,\Bigg{\langle}\bar{q}(0)\left[\overrightarrow{D}^{\alpha}-\overleftarrow{D}^{\alpha}\right]\gamma^{\beta}q(0)\Bigg{\rangle}^{\mu_{3},\,\mu_{4}}_{\overline{\mathrm{MS}}(\mu)}$
$\displaystyle={\textrm{finite
}}+i\,\frac{2}{3}\,N_{c}\,e_{q,j}^{2}\,\ln\frac{Q_{3}^{2}}{\mu^{2}}\lim_{q_{4}\rightarrow
0}\partial^{\nu_{3}}_{q_{3}}\partial^{\nu_{4}}_{q_{4}}\Bigg{\langle}\frac{\alpha}{4\pi}\left(F^{\mu\nu}F_{\mu\nu}g^{\alpha\beta}+d\,F^{\alpha\gamma}F_{\gamma}^{\beta}\right)\Bigg{\rangle}^{\mu_{3}\mu_{4}}\,.$
(80)

Writing
$\mathcal{O}^{\alpha\beta}_{j}=\bar{q}_{j}(0)\left[\overrightarrow{D}^{\alpha}-\overleftarrow{D}^{\alpha}\right]\gamma^{\beta}q_{j}(0)$
and denoting the bare operator by $\mathcal{O}^{\alpha\beta}_{(0),\,j}$, the
mixing is of the form

$\mathcal{O}^{\alpha\beta}_{(0),\,j}=\mathcal{O}^{\alpha\beta}_{j}(\mu)+Z^{j}_{DF}(\mu)\,\frac{\alpha}{4\pi}\,\left(F^{\mu\nu}F_{\mu\nu}g^{\alpha\beta}+d\,F^{\alpha\gamma}F_{\gamma}^{\beta}\right)+\ldots\,,$
(81)

where the renormalization factor $Z^{j}_{DF}(\mu)$ is

$Z^{j}_{DF}(\mu)=-i\,\frac{2}{3}\,N_{c}\,e_{q,j}^{2}\,\frac{\mu^{-2\epsilon}}{\hat{\epsilon}}\,.$
(82)

Repeating exactly the same calculation but with two initial-state gluons we
find

$\mathcal{O}^{\alpha\beta}_{(0),\,j}=\mathcal{O}^{\alpha\beta}_{j}(\mu)+Z^{j}_{DG}(\mu)\,\frac{\alpha_{s}}{4\pi}\,\left(G^{a,\mu\nu}G^{a}_{\mu\nu}\,g^{\alpha\beta}+d\,G^{\alpha\gamma}_{a}G_{\gamma}^{a,\,\beta}\right)+\ldots\,,$
(83)

with

$Z^{j}_{DG}(\mu)=-i\,\frac{1}{3}\,\frac{\mu^{-2\epsilon}}{\hat{\epsilon}}\,.$
(84)

These are needed to have finite Wilson coefficients associated to the field
strength contributions and quark-current matrix elements at the associated
leading order contributions in Eq. (33).

The renormalized version of Eq. (33) is then

$\displaystyle\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4}^{\mu_{4}}}=\sum_{j}\frac{e_{q,j}^{2}}{e^{2}}\,\lim_{q_{4}\rightarrow
0}\partial^{\nu_{4}}_{q_{4}}\,\Bigg{\langle}\bar{q}(0)\left[\Gamma^{\mu_{1}\mu_{2}}(-\hat{q})-\Gamma^{\mu_{2}\mu_{1}}(-\hat{q})\right]q(0)\Bigg{\rangle}^{j,\mu_{3},\mu_{4}}$
$\displaystyle+\sum_{j}\frac{ie_{q_{j}}^{2}}{e^{2}\hat{q}^{2}}\,\left(g^{\mu_{1}\delta}g^{\mu_{2}}_{\beta}+g^{\mu_{2}\delta}g^{\mu_{1}}_{\beta}-g^{\mu_{1}\mu_{2}}g^{\delta}_{\beta}\right)\left(g_{\alpha}^{\delta}-2\,\frac{\hat{q}^{\delta}\hat{q}_{\alpha}}{\hat{q}^{2}}\right)$
$\displaystyle\qquad\times\lim_{q_{4}\rightarrow
0}\partial_{\nu_{4}}^{q_{4}}\,\Bigg{\langle}\bar{q}(0)(\overrightarrow{D}^{\alpha}-\overleftarrow{D}^{\alpha})\gamma^{\beta}q(0)\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}_{\overline{\mathrm{MS}}(\mu)}$
$\displaystyle+\sum_{j}\frac{ie_{q_{j}}^{2}}{e^{2}\hat{q}^{2}}\,\left(g^{\mu_{1}\delta}g^{\mu_{2}}_{\beta}+g^{\mu_{2}\delta}g^{\mu_{1}}_{\beta}-g^{\mu_{1}\mu_{2}}g^{\delta}_{\beta}\right)\left(g_{\alpha}^{\delta}-2\,\frac{\hat{q}^{\delta}\hat{q}_{\alpha}}{\hat{q}^{2}}\right)$
$\displaystyle\qquad\times\lim_{q_{4}\rightarrow
0}\partial^{\nu_{4}}_{q_{4}}\,\Bigg{\langle}Z^{j}_{DF}(\mu)\,\frac{\alpha}{4\pi}\,\left(F^{\mu\nu}F_{\mu\nu}g^{\alpha\beta}+d\,F^{\alpha\gamma}F_{\gamma}^{\beta}\right)$
$\displaystyle\qquad\qquad\qquad+Z^{j}_{DG}(\mu)\,\frac{\alpha_{s}}{4\pi}\,\left(G^{\mu\nu}_{a}G_{\mu\nu}^{a}\,g^{\alpha\beta}+d\,G^{\alpha\gamma}_{a}G_{\gamma}^{a,\,\beta}\right)\Bigg{\rangle}^{j,\mu_{3}\mu_{4}}$
$\displaystyle+\sum_{j}\frac{e_{q_{j}}^{2}}{8e^{2}}\lim_{q_{4}\rightarrow
0}\left[\partial^{\nu_{4}}_{q_{4}}\,\Bigg{\langle}e^{2}e_{q_{j}}^{2}F_{\nu_{3}^{\prime}\mu_{3}^{\prime}}F_{\nu_{4}^{\prime}\mu_{4}^{\prime}}+\frac{1}{2N_{c}}\,g_{s}^{2}\,G^{a}_{\nu_{3}^{\prime}\mu_{3}^{\prime}}G^{a}_{\nu_{4}^{\prime}\mu_{4}^{\prime}}\Bigg{\rangle}^{j,\mu_{3}\mu_{4}}\right]\,\lim_{q_{3},q_{4}\rightarrow
0}\partial_{q_{3}}^{\nu_{3}^{\prime}}\,\partial_{q_{4}}^{\nu_{4}^{\prime}}\,\Pi_{\mathrm{ql},j}^{\mu_{1}\mu_{2}\mu_{3}^{\prime}\mu_{4}^{\prime}}\,,$
(85)

The last four lines together yield a renormalized Wilson coefficient. Within
this renormalization prescription the field strength contributions to the
$\hat{\Pi}_{i}$ are found to be zero, except for

$\displaystyle\underline{q_{3}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{4}$
$\displaystyle=-\frac{16}{3\pi^{2}\overline{Q}_{3}^{4}}\,,$
$\displaystyle\underline{q_{1}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\frac{8}{3\pi^{2}\overline{Q}_{1}^{4}}\left(1-4\log\frac{\overline{Q}_{1}}{2\mu}\right)\,,$
$\displaystyle\underline{q_{2}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\frac{8}{3\pi^{2}\overline{Q}_{2}^{4}}\left(1-4\log\frac{\overline{Q}_{2}}{2\mu}\right)\,.$
(86)

This can be compared to the non-renormalized case in Eqs. (4.2.2–4.2.2). The
remaining $\hat{\Pi}_{i}$ contributions are straightforward to obtain by
simply subtracting this contribution and the ones coming from $D=3$ to the
total $\hat{\Pi}_{i}$ in Eqs. (4.2.4–4.2.4). As a final remark, to the order
expanded the OPE in Eq. (85) satisfies gauge invariance in both $q_{1}$ and
$q_{2}$ as well.

##  5 Nonperturbative regime

In this section we extend the results of Sec. 4 beyond the
$Q_{3}^{2}\gg\Lambda_{\mathrm{QCD}}^{2}$ regime. We start by reproducing the
$D=3$ limit of Ref. [15]. For the power-suppressed $D=4$ extension we directly
work in the chiral limit, where we conveniently separate the tensor of Eq.
(16) in two contributions, a flavor singlet
$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}_{(1)}$ and a flavor octet
$\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}_{(8)}$, according to

$\displaystyle\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$
$\displaystyle=\sum_{j}\frac{ie_{q_{j}}^{2}}{e^{2}}\Bigg{\langle}e^{-i(q_{1}x_{1}+q_{2}x_{2})}\,\bar{q}(x_{1})\gamma^{\mu_{1}}q(x_{1})\;\bar{q}(x_{2})\gamma^{\mu_{2}}q(x_{2})\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}_{q_{4},\,x_{1},\,x_{2}}\,$
$\displaystyle=\frac{i\sum_{j}(e_{q_{j}}^{2}-\sum_{k}\frac{e_{q_{k}}^{2}}{3})}{e^{2}}\,\Bigg{\langle}e^{-i(q_{1}x_{1}+q_{2}x_{2})}\,\bar{q}(x_{1})\gamma^{\mu_{1}}q(x_{1})\;\bar{q}(x_{2})\gamma^{\mu_{2}}q(x_{2})\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}_{q_{4},\,x_{1},\,x_{2}}$
$\displaystyle+\frac{i\sum_{j}\sum_{k}\frac{e_{q_{k}}^{2}}{3}}{e^{2}}\Bigg{\langle}e^{-i(q_{1}x_{1}+q_{2}x_{2})}\,\bar{q}(x_{1})\gamma^{\mu_{1}}q(x_{1})\;\bar{q}(x_{2})\gamma^{\mu_{2}}q(x_{2})\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}_{q_{4},\,x_{1},\,x_{2}}$
$\displaystyle\equiv\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}_{(8)}+\Pi^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}_{(1)}\,.$
(87)

Obviously, we will have a much better predictive power for the octet case
because all gluon operator contributions (independent of $j$) cancel, but the
form of the obtained OPE is also valid for the singlet.

All $\hat{\Pi}_{i}$ and form factors are in this section given in units of
$N_{c}\sum_{j}e_{q_{j}}^{4}$.

###  5.1 Dimension $D=3$

Let us start by rewriting the dimension $D=3$ part of Eq. (33) as

$\displaystyle\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial q_{4,\,\mu_{4}}}$
$\displaystyle=-\sum_{j}\frac{e_{q,j}^{2}}{e^{2}}\,\frac{\hat{q}^{\alpha}}{\hat{q}^{2}}\,g^{\alpha_{1}\mu_{1}}g^{\alpha_{2}\mu_{2}}$
$\displaystyle\times\lim_{q_{4}\rightarrow
0}\partial_{q_{4}}^{\nu_{4}}\,\langle
0|\bar{q}_{j}(\gamma_{\alpha_{1}}\gamma_{\alpha}\gamma_{\alpha_{2}}-\gamma_{\alpha_{2}}\gamma_{\alpha}\gamma_{\alpha_{1}})q_{j}|\gamma^{\mu_{3}}(q_{3})\gamma^{\mu_{4}}(q_{4})\rangle\,.$
(88)

Using the relation555We work with the convention $\epsilon^{0123}=1$.

$\gamma_{\alpha_{1}}\gamma_{\alpha}\gamma_{\alpha_{2}}-\gamma_{\alpha_{2}}\gamma_{\alpha}\gamma_{\alpha_{1}}=-2i\,\epsilon_{\delta\alpha_{1}\alpha\alpha_{2}}\gamma^{\delta}\gamma_{5}\,,$
(89)

one finds

$\displaystyle\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial q_{4,\,\mu_{4}}}$
$\displaystyle=2\sum_{j}\frac{e_{q,j}^{2}}{e^{2}}\,\epsilon^{\mu_{1}\mu_{2}\alpha\delta}\,\frac{\hat{q}^{\alpha}}{\hat{q}^{2}}\lim_{q_{4}\rightarrow
0}i\,\partial_{q_{4}}^{\nu_{4}}\,\langle
0|\bar{q}_{j}\gamma^{\delta}\gamma_{5}q_{j}|\gamma^{\mu_{3}}(q_{3})\gamma^{\mu_{4}}(q_{4})\rangle.$
(90)

In the nonperturbative regime without an inherent ordering on $Q_{3}$ and
$\Lambda_{\textrm{QCD}}$, we must express the matrix element in a form factor
decomposition. The most general decomposition into Lorentz structures and
associated form factors is given by

$\displaystyle\sum_{j}\frac{e_{q,j}^{2}}{e^{2}}\lim_{q_{4}\rightarrow
0}\,i\,\partial_{q_{4}}^{\nu_{4}}\langle
0|\bar{q}_{j}\gamma^{\delta}\gamma_{5}q_{j}|\gamma^{\mu_{3}}(q_{3})\gamma^{\mu_{4}}(q_{4})\rangle$
$\displaystyle=\frac{1}{4\pi^{2}}\,q_{3}^{2}\,\left(\epsilon_{\mu_{3}\mu_{4}\nu_{4}\delta}\,\omega_{T}(q_{3}^{2})-\frac{1}{q_{3}^{2}}\,\epsilon_{q_{3}\mu_{4}\nu_{4}\delta}\,q_{3\mu_{3}}\,\omega_{T}(q_{3}^{2})+\frac{1}{q_{3}^{2}}\,\epsilon_{\mu_{3}\mu_{4}\nu_{4}q_{3}}\,q_{3\delta}\,\left[\omega_{L}(q_{3}^{2})-\omega_{T}(q_{3}^{2})\right]\right)\,.$
(91)

Here we introduced the compact notation $\epsilon_{\cdots
q\cdots}\equiv\epsilon_{\cdots\mu\cdots}q^{\mu}$, as well as the the
longitudinal and transversal form factors $\omega_{L,T}(q_{3}^{2})$ following
the notation in Ref. [15]. This yields the HLbL derivative tensor as

$\displaystyle\lim_{q_{4}\rightarrow 0}$
$\displaystyle\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}}{\partial
q_{4,\,\mu_{4}}}=\frac{1}{2\pi^{2}}\,\frac{q_{3}^{2}}{\hat{q}^{2}}\,\epsilon^{\mu_{1}\mu_{2}\hat{q}\delta}\Big{(}\epsilon_{\mu_{3}\mu_{4}\nu_{4}\delta}\,\omega_{T}(q_{3}^{2})-\frac{1}{q_{3}^{2}}\,\epsilon_{q_{3}\mu_{4}\nu_{4}\delta}\,q_{3\mu_{3}}\,\omega_{T}(q_{3}^{2})$
$\displaystyle\hskip
118.4024pt+\frac{1}{q_{3}^{2}}\,\epsilon_{\mu_{3}\mu_{4}\nu_{4}q_{3}}\,q_{3\delta}\,\left[\omega_{L}(q_{3}^{2})-\omega_{T}(q_{3}^{2})\right]\Big{)}\,.$
(92)

When $Q_{3}^{2}$ is large, i.e. in the perturbative regime, the partonic
description holds and

$\omega_{L}=2\,\omega_{T}=\sum_{j}\frac{2N_{c}e_{q_{j}}^{4}}{-Q_{3}^{2}}\,,$
(93)

which agrees with the $D=3$ result obtained in Sec. 4 and with Ref. [15].
Next, we can match to find the corresponding generalized $\hat{\Pi}_{i}$. In
the $q_{3}$ small region one finds at LO in the $\frac{1}{\overline{Q}_{3}}$
counting,

$\displaystyle\underline{q_{3}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\frac{2}{\pi^{2}\overline{Q}_{3}^{2}}\,\omega_{L}(q_{3}^{2})\,.$
(94)

This can be compared to Eq. (4.1). The partial contribution coming from $D=3$
at NLO in the $\frac{1}{\overline{Q}_{3}}$ counting, i.e. extending Eq. 4.2.3,
is given by

$\displaystyle\underline{q_{3}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{17}$
$\displaystyle=-\frac{4Q_{3}^{2}}{\pi^{2}(Q_{3}^{2}-\delta_{12}^{2})\overline{Q}_{3}^{4}}\,\omega_{T}(q_{3}^{2})\,,$
$\displaystyle\hat{\Pi}_{39}$
$\displaystyle=-\frac{4Q_{3}^{2}}{\pi^{2}(Q_{3}^{2}-\delta_{12}^{2})\overline{Q}_{3}^{4}}\,\omega_{T}(q_{3}^{2})\,.$
(95)

Through NLO in the $q_{1}$ small region we find

$\displaystyle\underline{q_{1}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{4}$
$\displaystyle=\frac{Q_{1}^{2}}{\pi^{2}(Q_{1}^{2}-\delta_{23}^{2})\overline{Q}_{1}^{2}}\,\omega_{T}(q_{1}^{2})\,,$
(96) $\displaystyle\hat{\Pi}_{7}$
$\displaystyle=-\frac{4Q_{1}^{2}\delta_{23}}{\pi^{2}(Q_{1}^{2}-\delta_{23}^{2})\overline{Q}_{1}^{3}}\,\omega_{T}(q_{1}^{2})\,,$
(97) $\displaystyle\hat{\Pi}_{39}$
$\displaystyle=-\frac{4Q_{1}^{2}}{\pi^{2}(Q_{1}^{2}-\delta_{23}^{2})\overline{Q}_{1}^{4}}\,\omega_{T}(q_{1}^{2})\,,$
(98) $\displaystyle\hat{\Pi}_{54}$
$\displaystyle=-\frac{4Q_{1}^{2}}{\pi^{2}(Q_{1}^{2}-\delta_{23}^{2})\overline{Q}_{1}^{4}}\,\omega_{T}(q_{1}^{2})\,.$
(99)

Finally in the $q_{2}$ small region the result is

$\displaystyle\underline{q_{2}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{4}$
$\displaystyle=\frac{Q_{2}^{2}}{\pi^{2}(Q_{2}^{2}-\delta_{31}^{2})\overline{Q}_{2}^{2}}\,\omega_{T}(q_{2}^{2})\,,$
(100) $\displaystyle\hat{\Pi}_{7}$
$\displaystyle=\frac{4Q_{2}^{2}}{\pi^{2}(Q_{2}^{2}-\delta_{31}^{2})\overline{Q}_{2}^{4}}\,\omega_{T}(q_{2}^{2})\,,$
(101) $\displaystyle\hat{\Pi}_{39}$
$\displaystyle=-\frac{4Q_{2}^{2}}{\pi^{2}(Q_{2}^{2}-\delta_{31}^{2})\overline{Q}_{2}^{4}}\,\omega_{T}(q_{2}^{2})\,,$
(102) $\displaystyle\hat{\Pi}_{54}$
$\displaystyle=\frac{4Q_{2}^{2}}{\pi^{2}(Q_{2}^{2}-\delta_{31}^{2})\overline{Q}_{2}^{4}}\,\omega_{T}(q_{2}^{2})\,.$
(103)

From the collinear divergences $\delta_{ij}\rightarrow Q_{k}$ in these
expressions it is once again transparent how keeping only the $D=3$ at NLO in
$\frac{1}{\overline{Q}_{k}}$ is not consistent: the $\hat{\Pi}_{i}$ are by
construction free from kinematic singularities while these partial results are
not. It can easily be checked that by saturating the form factors to their
short-distance limits of Eq. (93), one recovers the corresponding
$\hat{\Pi}_{i}$ in Eqs. (4.2.3–4.2.3).

Even if it is not the focus of this work and there are many of them about the
topic [44, 45, 46, 47, 34], let us note how in the chiral limit further
simplifications can be made. First, the third and fourth photon must emerge
from the same quark line (and thus from the same flavor), otherwise one would
have a vanishing flavor trace, $\mathrm{Tr}(e_{q})=0$. We then have a
disconnected and a connected topology. For the same reason, the octet piece
$\Pi_{(8)}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$, does not receive any contribution
from disconnected topologies.

The octet form factors, $\omega_{L,(8)}$ and $\omega_{T,(8)}$, are the ones
satisfying nonrenormalization theorems [45, 15, 48]. Concretely, the
longitudinal octet value given by

$\displaystyle\omega_{L,(8)}=\frac{\sum_{j}e_{q,j}^{4}-\frac{1}{3}\sum_{j,k}e_{q,j}^{2}e_{q,k}^{2}}{\sum_{j}e_{q,j}^{4}}=\frac{1}{3}\omega_{L}\,,$
(104)

from Eq. 93 is exact beyond perturbation theory, i.e. for all $Q_{3}^{2}$
values, and $\omega_{T,(8)}$ does not receive any perturbative correction.
However, the validity of Eq. (93) cannot be extended beyond perturbation
theory, i.e. as $Q_{3}^{2}\rightarrow\Lambda_{\mathrm{QCD}}^{2}$. While much
less can be said about the singlet, one naively expects its contribution to be
suppressed with respect to the octet because there are no zero mass modes in
the singlet axial channel.

###  5.2 Dimension $D=4$

As explained above, for the suppressed $D=4$ contribution we directly work in
the chiral limit. In Table 1 we summarise which contributions appear at
dimension $D=4$. This follows since the disconnected topology vanishes for
$\Pi_{(8)}^{\mu_{1}\mu_{2}\mu_{3}\mu_{4}}$ and that the gluon operator and the
photon one (at this order in $\alpha_{s}(\overline{Q}_{3})$) respectively are
disconnected and connected.

Operator | $\Pi_{(1)}$ | $\Pi_{(8)}$
---|---|---
Connected | Disconnected | Connected | Disconnected
Photon | ✔ | ✘ | ✔ | ✘
Gluon | ✘ | ✔ | ✘ | ✘
Quark | ✔ | ✔ | ✔ | ✘
Table 1: Summary of studied contributions entering at dimension $D=4$.

####  5.2.1 Octet case

For the octet case the gluon operator contribution trivially drops due to its
flavorless nature. The photon matrix element does not receive nonperturbative
correction at the studied order in $\alpha$ and it is therefore identical to
the contribution studied in Sec. 4, except for a trivial charge factor. The
associated NLO contribution to the OPE then becomes

$\displaystyle\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}_{(8)}}{\partial
q_{4,\,\mu_{4}}}=\frac{1}{3}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}_{FF,\mu}}{\partial
q_{4,\,\mu_{4}}}$
$\displaystyle+\frac{i\sum_{j}\left(e_{q_{j}}^{2}-\sum_{k}\frac{e_{q,k}^{2}}{3}\right)}{e^{2}\hat{q}^{2}}\,\left(g^{\mu_{1}}_{\delta}g^{\mu_{2}}_{\beta}+g^{\mu_{2}}_{\delta}g^{\mu_{1}}_{\beta}-g^{\mu_{1}\mu_{2}}g_{\delta\beta}\right)\left(g_{\alpha}^{\delta}-2\,\frac{\hat{q}^{\delta}\hat{q}_{\alpha}}{\hat{q}^{2}}\right)$
$\displaystyle\times\lim_{q_{4}\rightarrow
0}\partial^{\nu_{4}}_{q_{4}}\,\Bigg{\langle}\bar{q}(0)\left[\overrightarrow{D}^{\alpha}-\overleftarrow{D}^{\alpha}\right]\gamma^{\beta}q(0)\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}_{\overline{\mathrm{MS}}(\mu)}\,.$

where $\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}_{FF,\mu}$ simply refers to the
Wilson-renormalized photon operator contribution in Eq. (85), whose
contribution to the $\hat{\Pi}_{i}$ is still given by Eq. (4.2.5). The
remaining matrix element with 5 Lorentz indices
$\alpha\beta\mu_{3}\mu_{4}\nu_{4}$ can be decomposed into form factors. There
are in total 26 possible Lorentz structures from combinations of $q_{3}^{\mu}$
and $g^{\nu\rho}$ (the Levi-Civita tensor does not appear because of parity).
However, we know that the final matrix element has to be anti-symmetric in
$\mu_{4}\leftrightarrow\nu_{4}$ and be gauge invariant with respect to
$q_{3}^{\mu_{3}}$. Applying these conditions yields a set of 6 independent
form factors $\omega_{(8)}^{D,i}(q_{3}^{2})$ so that the matrix element can be
written

$\frac{i\sum_{j}\left(e_{q_{j}}^{2}-\sum_{k}\frac{e_{q,k}^{2}}{3}\right)}{e^{2}\hat{q}^{2}}\lim_{q_{4}\rightarrow
0}\,\partial^{\nu_{4}}_{q_{4}}\,\Bigg{\langle}\bar{q}(0)\left[\overrightarrow{D}^{\alpha}-\overleftarrow{D}^{\alpha}\right]\gamma^{\beta}q(0)\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}_{\overline{\textrm{MS}}(\mu)}=\sum_{i=1}^{6}\omega_{(8)}^{D,i}\,L_{i}^{\alpha\beta\mu_{3}\mu_{4}\nu_{4}}\,.$
(105)

The Lorentz structures are

$\displaystyle L^{\alpha\beta\mu_{3}\mu_{4}\nu_{4}}_{1}$
$\displaystyle=g^{\mu_{3}\mu_{4}}q_{3}^{\nu_{4}}g^{\alpha\beta}-(\mu_{4}\leftrightarrow\nu_{4})\,,$
$\displaystyle L^{\alpha\beta\mu_{3}\mu_{4}\nu_{4}}_{2}$
$\displaystyle=g^{\beta\mu_{4}}q_{3}^{\nu_{4}}\left(g^{\alpha\mu_{3}}-\frac{q_{3}^{\alpha}q_{3}^{\mu_{3}}}{q_{3}^{2}}\right)-(\mu_{4}\leftrightarrow\nu_{4})\,,$
$\displaystyle L^{\alpha\beta\mu_{3}\mu_{4}\nu_{4}}_{3}$
$\displaystyle=L^{\beta\alpha\mu_{3}\mu_{4}\nu_{4}}_{2}\,,$ $\displaystyle
L^{\alpha\beta\mu_{3}\mu_{4}\nu_{4}}_{4}$
$\displaystyle=g^{\alpha\mu_{4}}q_{3}^{\beta}\left(g^{\mu_{3}\nu_{4}}-\frac{q_{3}^{\mu_{3}}q_{3}^{\nu_{4}}}{q_{3}^{2}}\right)-(\mu_{4}\leftrightarrow\nu_{4})\,,$
$\displaystyle L^{\alpha\beta\mu_{3}\mu_{4}\nu_{4}}_{5}$
$\displaystyle=L^{\beta\alpha\mu_{3}\mu_{4}\nu_{4}}_{4}\,,$ $\displaystyle
L^{\alpha\beta\mu_{3}\mu_{4}\nu_{4}}_{6}$
$\displaystyle=g^{\mu_{3}\mu_{4}}q_{3}^{\nu_{4}}q^{\alpha}_{3}q^{\beta}_{3}-(\mu_{4}\leftrightarrow\nu_{4})\,.$
(106)

Projecting onto the $\hat{\Pi}_{i}$ and requiring that they are free of
kinematic singularities, we obtain nontrivial relations among nonperturbative
form factors. The conditions are666Let us note how the
$\omega_{(8)}^{D,4}=\omega_{(8)}^{D,5}$ and the
$\omega_{(8)}^{D,3}-\omega_{(8)}^{D,2}=\frac{Q_{i}^{2}}{4\pi^{2}}\omega_{T,(8)}$
relations are constraints on the antisymmetric part of the quark operator,
which at the studied order does not mix with the photon one.

$\displaystyle\omega_{(8)}^{D,2}$
$\displaystyle=-2\,\omega_{(8)}^{D,1}+\omega_{(8)}^{D,5}-\frac{\omega_{(8)}^{D,6}}{2}-\frac{\omega_{T,(8)}Q_{i}^{2}}{8\pi^{2}}\,,$
$\displaystyle\omega_{(8)}^{D,3}$
$\displaystyle=-2\,\omega_{(8)}^{D,1}+\omega_{(8)}^{D,5}-\frac{\omega_{(8)}^{D,6}}{2}+\frac{\omega_{T,(8)}Q_{i}^{2}}{8\pi^{2}}\,,$
$\displaystyle\omega_{(8)}^{D,4}$ $\displaystyle=\omega_{(8)}^{D,5}\,.$ (107)

Using these relations, we find the $\hat{\Pi}_{i}$ to be

$\displaystyle\underline{q_{3}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{4}}\right)\,,$
(108) $\displaystyle\hat{\Pi}_{4}$
$\displaystyle=-\frac{64\left(\omega_{(8)}^{D,1}-\omega_{(8)}^{D,5}\right)}{\overline{Q}_{3}^{4}}\,,$
(109) $\displaystyle\hat{\Pi}_{7}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{6}}\right)\,,$
(110) $\displaystyle\hat{\Pi}_{17}$
$\displaystyle=\frac{32\left(-2\omega_{(8)}^{D,1}+2\omega_{(8)}^{D,5}+\frac{\omega_{(8)}^{D,6}}{2}-\frac{\omega_{T}Q_{3}^{2}}{8\pi^{2}}\right)}{\overline{Q}_{3}^{4}Q_{3}^{2}}\,,\,,$
(111) $\displaystyle\hat{\Pi}_{39}$
$\displaystyle=-\frac{32\left(-2\omega_{(8)}^{D,1}+2\omega_{(8)}^{D,5}+\frac{\omega_{(8)}^{D,6}}{2}+\frac{\omega_{T}Q_{3}^{2}}{8\pi^{2}}\right)}{\overline{Q}_{3}^{4}Q_{3}^{2}}\,,$
(112) $\displaystyle\hat{\Pi}_{54}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{5}}\right)\,.$
(113) $\displaystyle\underline{q_{1}\textrm{ small:}}$
$\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\frac{32\,\omega_{(8)}^{D,5}}{\overline{Q}_{1}^{4}}\,,$ (114)
$\displaystyle\hat{\Pi}_{4}$
$\displaystyle=-\frac{16\,\omega_{(8)}^{D,1}-16\,\omega_{(8)}^{D,5}+4\,\omega_{(8)}^{D,6}-\frac{\omega_{T,(8)}Q_{1}^{2}}{\pi^{2}}}{Q_{1}^{2}\overline{Q}_{1}^{2}}\,,$
(115) $\displaystyle\hat{\Pi}_{7}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{4}}\right)\,,$
(116) $\displaystyle\hat{\Pi}_{17}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{5}}\right)\,,$
(117) $\displaystyle\hat{\Pi}_{39}$
$\displaystyle=\frac{4\left(16\,\omega_{(8)}^{D,1}-16\,\omega_{(8)}^{D,5}+4\omega_{(8)}^{D,6}-\frac{\omega_{T,(8)}Q_{1}^{2}}{\pi^{2}}\right)}{Q_{1}^{2}\overline{Q}_{1}^{2}}\,,$
(118) $\displaystyle\hat{\Pi}_{54}$
$\displaystyle=-\frac{4\left(16\,\omega_{(8)}^{D,1}-16\,\omega_{(8)}^{D,5}+4\,\omega_{(8)}^{D,6}+\frac{\omega_{T,(8)}Q_{1}^{2}}{\pi^{2}}\right)}{Q_{1}^{2}\overline{Q}_{1}^{2}}\,.$
(119) $\displaystyle\underline{q_{2}\textrm{ small:}}$
$\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\frac{32}{\overline{Q}_{2}^{4}}\,\omega_{(8)}^{D,5}\,,$ (120)
$\displaystyle\hat{\Pi}_{4}$
$\displaystyle=-\frac{16\,\omega_{(8)}^{D,1}-16\,\omega_{(8)}^{D,5}+4\,\omega_{(8)}^{D,6}-\frac{\omega_{T,(8)}Q_{2}^{2}}{\pi^{2}}}{Q_{2}^{2}\overline{Q}_{2}^{2}}\,,$
(121) $\displaystyle\hat{\Pi}_{7}$
$\displaystyle=\frac{16\,\omega_{(8)}^{D,1}-16\,\omega_{(8)}^{D,5}+4\,\omega_{(8)}^{D,6}-\frac{\omega_{T,(8)}Q_{2}^{2}}{\pi^{2}}}{Q_{2}^{2}\overline{Q}_{2}^{2}}\,,$
(122) $\displaystyle\hat{\Pi}_{17}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{5}}\right)\,,$
(123) $\displaystyle\hat{\Pi}_{39}$
$\displaystyle=\frac{16\,\omega_{(8)}^{D,1}-16\,\omega_{(8)}^{D,5}+4\,\omega_{(8)}^{D,6}-\frac{\omega_{T,(8)}Q_{2}^{2}}{\pi^{2}}}{Q_{2}^{2}\overline{Q}_{2}^{2}}\,,$
(124) $\displaystyle\hat{\Pi}_{54}$
$\displaystyle=\frac{16\,\omega_{(8)}^{D,1}-16\,\omega_{(8)}^{D,5}+4\,\omega_{(8)}^{D,6}+\frac{\omega_{T,(8)}Q_{2}^{2}}{\pi^{2}}}{Q_{2}^{2}\overline{Q}_{2}^{2}}\,.$
(125)

In fact, me may here take the perturbative limits
$Q_{i}^{2}\gg\Lambda_{\textrm{QCD}}^{2}$ and match the three appearing form
factors $\omega_{(8)}^{D,1}$, $\omega_{(8)}^{D,5}$ and $\omega_{(8)}^{D,6}$.
Using Eq. (93), the result is

$\displaystyle\underline{q_{i}\textrm{ small:}}$
$\displaystyle\omega_{(8)}^{D,1}$
$\displaystyle=\frac{-13+12\log\frac{Q_{i}}{\mu}}{108\,\pi^{2}}\,,$
$\displaystyle\omega_{(8)}^{D,5}$
$\displaystyle=\frac{-13+12\log\frac{Q_{i}}{\mu}}{108\,\pi^{2}}\,,$
$\displaystyle\omega_{(8)}^{D,6}$ $\displaystyle=\frac{7}{36\,\pi^{2}}\,.$
(126)

The full NLO octet result at the nonperturbative level becomes then
parameterized in terms of $3$ form factors depending on a single energy scale,
whose short-distance dependence is known and with a clear definition, given by
Eq. (105). Complementing this model-independent result with nonperturbative
methods to extrapolate them to lower energies, such as resonance models,
should allow for a more reliable estimate of the subleading contributions to
$g-2$ coming from this particularly challenging domain.

####  5.2.2 Singlet case

More complicated is the singlet case, since the gluon operator does not drop
out. In this case the NLO piece takes the form

$\displaystyle\lim_{q_{4}\rightarrow
0}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}_{(1)}}{\partial
q_{4,\,\mu_{4}}}=\frac{2}{3}\frac{\partial\Pi^{\mu_{1}\mu_{2}\mu_{3}\nu_{4}}_{FF,\mu}}{\partial
q_{4}^{\mu_{4}}}$
$\displaystyle+i\,\frac{\sum_{j}\sum_{k}\frac{e_{k}^{2}}{3}}{e^{2}\hat{q}^{2}}\,\left(g^{\mu_{1}}_{\delta}g^{\mu_{2}}_{\beta}+g^{\mu_{2}}_{\delta}g^{\mu_{1}}_{\beta}-g^{\mu_{1}\mu_{2}}g_{\delta\beta}\right)\left(g_{\alpha}^{\delta}-2\,\frac{\hat{q}^{\delta}\hat{q}_{\alpha}}{\hat{q}^{2}}\right)$
$\displaystyle\qquad\times\lim_{q_{4}\rightarrow
0}\partial^{\nu_{4}}_{q_{4}}\,\Bigg{\langle}\bar{q}(0)\left[\overrightarrow{D}^{\alpha}-\overleftarrow{D}^{\alpha}\right]\gamma^{\beta}q(0)\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}_{\overline{\textrm{MS}}(\mu)}$
$\displaystyle+\frac{i\,\sum_{j}\sum_{k}\frac{e_{q_{k}}^{2}}{3}}{e^{2}\hat{q}^{2}}\left(g^{\mu_{1}}_{\delta}g^{\mu_{2}}_{\beta}+g^{\mu_{2}}_{\delta}g^{\mu_{1}}_{\beta}-g^{\mu_{1}\mu_{2}}g_{\delta\beta}\right)\left(g_{\alpha}^{\delta}-2\,\frac{\hat{q}^{\delta}\hat{q}_{\alpha}}{\hat{q}^{2}}\right)$
$\displaystyle\qquad\times\lim_{q_{4}\rightarrow
0}\partial^{\nu_{4}}_{q_{4}}\,\left[Z^{j}_{DG}(\mu)\,\frac{\alpha_{s}}{4\pi}\,\left(G^{\mu\nu}_{a}G_{\mu\nu}^{a}\,g^{\alpha\beta}+d\,G^{\alpha\gamma}_{a}G^{a,\,\beta}_{\gamma}\right)\right]$
$\displaystyle+\frac{\sum_{j}\sum_{k}\frac{e_{q_{k}}^{2}}{3}}{8e^{2}}\,\lim_{q_{4}\rightarrow
0}\partial^{\nu_{4}}_{q_{4}}\,\left[\Bigg{\langle}\frac{1}{2N_{c}}\,g_{s}^{2}\,G^{a}_{\nu_{3}^{\prime}\mu_{3}^{\prime}}G^{a}_{\nu_{4}^{\prime}\mu_{4}^{\prime}}\Bigg{\rangle}^{j,\mu_{3}\mu_{4}}\right]\,\times\lim_{q_{3},q_{4}\rightarrow
0}\partial_{q_{3}}^{\nu_{3}^{\prime}}\partial_{q_{4}}^{\nu_{4}^{\prime}}\Pi_{\mathrm{ql,j}}^{\mu_{1}\mu_{2}\mu_{3}^{\prime}\mu_{4}^{\prime}}\,.$
(127)

We may decompose the renormalized quark-current matrix element in precisely
the same way, only now with 6 singlet form factors according to

$\frac{i\sum_{j}\sum_{k}\frac{e_{q_{k}}^{2}}{3}}{e^{2}\hat{q}^{2}}\lim_{q_{4}\rightarrow
0}\partial_{\nu_{4}}^{q_{4}}\Bigg{\langle}\bar{q}(0)\left[\overrightarrow{D}^{\alpha}-\overleftarrow{D}^{\alpha}\right]\gamma^{\beta}q(0)\Bigg{\rangle}^{j,\mu_{3},\,\mu_{4}}_{\overline{\textrm{MS}}(\mu)}=\sum_{i=1}^{6}\omega_{(1)}^{D,i}\,L_{i}^{\alpha\beta\mu_{3}\mu_{4}\nu_{4}}\,.$
(128)

The decomposition into $\hat{\Pi}_{i}$ is there identical to the octet case,
but their perturbative short-distance quark loop limit is twice the octet
case. This means that the for the singlet form factors $\omega_{(1)}^{D,1}$,
$\omega_{(1)}^{D,5}$ and $\omega_{(1)}^{D,6}$ , the perturbative limit is
instead of Eq. (5.2.1) now

$\displaystyle\underline{q_{i}\textrm{ small:}}$
$\displaystyle\omega_{(1)}^{D,1}$
$\displaystyle=\frac{-13+12\log\frac{Q_{i}}{\mu}}{54\,\pi^{2}}\,,$
$\displaystyle\omega_{(1)}^{D,5}$
$\displaystyle=\frac{-13+12\log\frac{Q_{i}}{\mu}}{54\,\pi^{2}}\,,$
$\displaystyle\omega_{(1)}^{D,6}$ $\displaystyle=\frac{7}{18\,\pi^{2}}\,.$
(129)

Of course, at higher orders and at the nonperturbative level
$\omega_{(1)}^{D,i}\neq 2\,\omega_{(8)}^{D,i}$.

The main qualitative difference for the singlet piece arises from the gluon
operator contributions. There is no $\alpha_{s}^{2}$ suppression when
$Q_{3}\rightarrow\Lambda_{\mathrm{QCD}}$.777One explicit example is the low-
energy theorem that relates $\langle
0|\alpha_{s}\tilde{G}G|\gamma\gamma\rangle$ to $\langle
0|\alpha\tilde{F}F|\gamma\gamma\rangle$ [49]. Once again the associated matrix
element needs to be form factor decomposed. There are in total 7 independent
form factors $\omega^{G}_{i}$ with associated Lorentz structures
$L_{i,G}^{\mu_{3}^{\prime}\mu_{4}^{\prime}\nu_{3}^{\prime}\nu_{4}^{\prime}\mu_{3}\mu_{4}\nu_{4}}$,
so that

$\displaystyle\sum_{j}\sum_{k}\frac{e_{q_{k}}^{2}}{3}\lim_{q_{4}\rightarrow
0}\partial^{\nu_{4}}_{q_{4}}\Bigg{\langle}\frac{1}{2e^{2}N_{c}}\,g_{s}^{2}\,G_{a}^{\nu_{3}^{\prime}\mu_{3}^{\prime}}G_{a}^{\nu_{4}^{\prime}\mu_{4}^{\prime}}\Bigg{\rangle}^{j,\mu_{3}\mu_{4}}=\sum_{i=1}^{7}\omega^{G}_{i}L_{i,G}^{\mu_{3}^{\prime}\mu_{4}^{\prime}\nu_{3}^{\prime}\nu_{4}^{\prime}\mu_{3}\mu_{4}\nu_{4}}\,.$
(130)

That only $7$ form factors appears for a matrix element with 7 Lorentz indices
follows from substantial symmetry requirements. The ones needed are

  * •

Anti-symmetry under $\mu_{3}^{\prime}\leftrightarrow\nu_{3}^{\prime}$
exchange.

  * •

Anti-symmetry under $\mu_{4}^{\prime}\leftrightarrow\nu_{4}^{\prime}$
exchange.

  * •

Anti-symmetry under $\mu_{4}\leftrightarrow\nu_{4}$ exchange.

  * •

Symmetry under simultaneous exchange of
$\mu_{3}^{\prime}\leftrightarrow\nu_{3}^{\prime}$ and
$\mu_{4}^{\prime}\leftrightarrow\nu_{4}^{\prime}$.

  * •

Gauge invariance condition under $\mu_{3}$.

The expressions for the Lorentz structures are rather lengthy and have been
relegated to the supplementary file singletbasis.txt. The decomposition has
been arranged in such a way that the analogous matrix element for the photon
case saturates to the first form factor, $\omega^{G}_{1}$, which additionally,
for the gluon case, can be linked to $\langle
0|\alpha_{s}\tilde{G}G|\gamma\gamma\rangle$, whose low-energy limit is known
[49]. Requiring that there are no kinematic singularities for the $\hat{\Pi}$
leads to

$\omega^{G}_{4}=0\,.$ (131)

Employing this condition, we find the $\hat{\Pi}_{i}$ in the three regions to
be

$\displaystyle\underline{q_{3}\textrm{ small:}}$ $\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{4}}\right)\,,$
(132) $\displaystyle\hat{\Pi}_{4}$
$\displaystyle=-\frac{16}{3\pi^{2}\overline{Q}_{3}^{4}}\Big{[}\omega_{1}^{G}-3\omega_{2}^{G}+\omega_{3}^{G}-\omega_{5}^{G}-2\omega_{6}^{G}-2\omega_{7}^{G}$
$\displaystyle\qquad-2\log\frac{\overline{Q}_{3}}{2\mu}\left(\omega_{5}^{G}+2\omega_{6}^{G}+2\omega_{7}^{G}\right)\Big{]}$
$\displaystyle+\frac{32}{3\pi^{2}\overline{Q}_{3}^{4}}\left(\omega_{5}^{G}+2\omega_{6}^{G}+2\omega_{7}^{G}\right)\left(1-\frac{\delta_{12}^{2}}{Q_{3}^{2}}\right)\,,$
(133) $\displaystyle\hat{\Pi}_{7}$
$\displaystyle=\frac{64\,\delta_{12}\,\left(\omega_{6}^{G}+2\omega_{7}^{G}\right)}{3\pi^{2}Q_{3}^{2}\overline{Q}_{3}^{5}}\,,$
(134) $\displaystyle\hat{\Pi}_{17}$
$\displaystyle=\frac{8}{3\pi^{2}Q_{3}^{2}\overline{Q}_{3}^{4}}\Bigg{[}\left(-3+8\log\frac{\overline{Q}_{3}}{2\mu}\right)\,\omega_{5}^{G}-3\,\omega_{6}^{G}+12\,\omega_{6}^{G}\,\log\frac{\overline{Q}_{3}}{2\mu}$
$\displaystyle\qquad-3\,\omega_{7}^{G}+8\,\omega_{7}^{G}\,\log\frac{\overline{Q}_{3}}{2\mu}\Bigg{]}\,,$
(135) $\displaystyle\hat{\Pi}_{39}$
$\displaystyle=-\frac{8}{3\pi^{2}Q_{3}^{2}\overline{Q}_{3}^{4}}\Bigg{[}\left(-1+4\log\frac{\overline{Q}_{3}}{2\mu}\right)\omega_{6}^{G}+\left(-3+8\log\frac{\overline{Q}_{3}}{2\mu}\right)\omega_{7}^{G}\Bigg{]}\,,$
(136) $\displaystyle\hat{\Pi}_{54}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{3}^{5}}\right)\,\,,$
(137) $\displaystyle\underline{q_{1}\textrm{ small:}}$
$\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\frac{8}{3\overline{Q}_{1}^{4}\pi^{2}}\left[\left(1-4\log\frac{\overline{Q}_{1}}{2\mu}\right)\left(\omega_{1}^{G}-3\omega_{2}^{G}\right)+\left(3-8\log\frac{\overline{Q}_{1}}{2\mu}\right)\omega_{3}^{G}\right]\,,$
(138) $\displaystyle\hat{\Pi}_{4}$
$\displaystyle=\frac{2}{3\overline{Q}_{1}^{2}Q_{1}^{2}\pi^{2}}\left[-\left(1-4\log\frac{\overline{Q}_{1}}{2\mu}\right)\omega_{6}^{G}-\left(3-8\log\frac{\overline{Q}_{1}}{2\mu}\right)\omega_{7}^{G}\right]\,,$
(139) $\displaystyle\hat{\Pi}_{7}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{4}}\right)\,,$
(140) $\displaystyle\hat{\Pi}_{17}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{1}^{5}}\right)\,,$
(141) $\displaystyle\hat{\Pi}_{39}$
$\displaystyle=-\frac{8}{3\overline{Q}_{1}^{4}Q_{1}^{2}\pi^{2}}\left[-\left(1-4\log\frac{\overline{Q}_{1}}{2\mu}\right)\omega_{6}^{G}-\left(3-8\log\frac{\overline{Q}_{1}}{2\mu}\right)\omega_{7}^{G}\right]\,,$
(142) $\displaystyle\hat{\Pi}_{54}$
$\displaystyle=\frac{8}{3\overline{Q}_{1}^{4}Q_{1}^{2}\pi^{2}}\left[-\left(1-4\log\frac{\overline{Q}_{1}}{2\mu}\right)\omega_{6}^{G}-\left(3-8\log\frac{\overline{Q}_{1}}{2\mu}\right)\omega_{7}^{G}\right]\,,$
(143) $\displaystyle\underline{q_{2}\textrm{ small:}}$
$\displaystyle\hat{\Pi}_{1}$
$\displaystyle=\frac{8}{3\overline{Q}_{2}^{4}\pi^{2}}\left[\left(1-4\log\frac{\overline{Q}_{2}}{2\mu}\right)\left(\omega_{1}^{G}-3\omega_{2}^{G}\right)+\left(3-8\log\frac{\overline{Q}_{2}}{2\mu}\right)\omega_{3}^{G}\right]\,,$
(144) $\displaystyle\hat{\Pi}_{4}$
$\displaystyle=\frac{2}{3\overline{Q}_{2}^{2}Q_{2}^{2}\pi^{2}}\left[-\left(1-4\log\frac{\overline{Q}_{2}}{2\mu}\right)\omega_{6}^{G}-\left(3-8\log\frac{\overline{Q}_{2}}{2\mu}\right)\omega_{7}^{G}\right]\,,$
(145) $\displaystyle\hat{\Pi}_{7}$
$\displaystyle=\frac{8}{3\overline{Q}_{2}^{4}Q_{2}^{2}\pi^{2}}\left[-\left(1-4\log\frac{\overline{Q}_{2}}{2\mu}\right)\omega_{6}^{G}-\left(3-8\log\frac{\overline{Q}_{2}}{2\mu}\right)\omega_{7}^{G}\right]\,,$
(146) $\displaystyle\hat{\Pi}_{17}$
$\displaystyle=\mathcal{O}\left(\frac{1}{\overline{Q}_{2}^{5}}\right)\,,$
(147) $\displaystyle\hat{\Pi}_{39}$
$\displaystyle=-\frac{8}{3\overline{Q}_{2}^{4}Q_{2}^{2}\pi^{2}}\left[-\left(1-4\log\frac{\overline{Q}_{2}}{2\mu}\right)\omega_{6}^{G}-\left(3-8\log\frac{\overline{Q}_{2}}{2\mu}\right)\omega_{7}^{G}\right]\,,$
(148) $\displaystyle\hat{\Pi}_{54}$
$\displaystyle=-\frac{8}{3\overline{Q}_{2}^{4}Q_{2}^{2}\pi^{2}}\left[-\left(1-4\log\frac{\overline{Q}_{2}}{2\mu})\omega_{6}^{G}-(3-8\log\frac{\overline{Q}_{2}}{2\mu}\right)\omega_{7}^{G}\right]\,.$
(149)

Let us note how the $\mu$-dependence on the Wilson coefficients multiplying
the gluon matrix elements cancel with the corresponding one of the $D=4$
matrix element involving quark fields. Notice how this is not in contradiction
with the fact that no logarithms appeared in the explicit
$Q_{3}\gg\Lambda_{\mathrm{QCD}}$ matrix elements computed in Sec. 4, since in
that regime the gluonic operator contributions become suppressed by
$\alpha^{2}_{s}(Q_{3}^{2})$ and then they are beyond the computed order.

##  6 Conclusions

In this paper we have studied short-distance constraints for the HLbL in the
so-called Melnikov-Vainshtein limit, where two of the photon virtualities are
greater than the third, or, $Q_{i}^{2},Q_{j}^{2}\gg
Q_{k}^{2},\Lambda_{\textrm{ QCD}}^{2}$. The constraints were derived using a
systematic and renormalized operator product expansion, through dimension
$D=4$ at leading order in $\alpha_{s}$, thus extending the results of Ref.
[15] and shedding light on the validity of the OPE. We have further shown that
in the limiting case $Q_{k}^{2}\gg\Lambda_{\mathrm{QCD}}^{2}$, our previous
results in Refs. [1, 2, 3] are reproduced. This completes our fundamental
understanding on the validity domains of the quark-gluon picture of the HLbL
tensor.

In the nonperturbative limit without an inherent ordering on $Q_{k}^{2}$ and
$\Lambda_{\mathrm{QCD}}^{2}$, we decomposed all matrix elements in general
form factor decompositions. At dimension $D=3$ we reproduce the known
longitudinal and transversal form factors of Ref. [15]. At dimension $D=4$, we
find, to the best of our knowledge, novel relations among Green functions of
different dimensionality by using that the $\hat{\Pi}_{i}$ by construction are
free from kinematic singularities. Once this is done, we find that in the
octet part all our perturbative ignorance becomes parameterized in terms of
three form factors depending on a single energy scale, whose short-distance
dependence is known and with a clear definition in terms of Green functions.
This opens a clear path to a more reliable determination of the corresponding
(subleading) $g-2$ contribution. For the singlet case six additional gluon
matrix elements, suppressed by $\alpha_{s}^{2}$ in the perturbative regime,
emerge, which may limit the associated predictive power. We leave the
nonperturbative completion of the form factors, through model estimates, as
well as its impact on $g-2$, for future work. Let us note that the
phenomenological impact of the corresponding constraints on
$a_{\mu}^{\textrm{HLbL}}$ has already been studied for the $D=3$ cases in a
number of papers [31, 19, 30, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41], also by
implementing them in models.

Let us finally point out that we are currently finalizing work on gluonic
corrections to our operator product expansion, which will be the topic of a
separate paper. For the dimension $D=3$ contribution, we confirm the statement
of Ref. [33] that it comes with a simple overall coefficient
$-\alpha_{s}/\pi$. At dimension $D=4$ we have also reproduced the corner
expansions in our previous paper Ref. [3].

## Acknowledgements

J. B. is supported by the Swedish Research Council grants contract numbers
2016-05996 and 2019-03779. N. H.-T. is funded in part by the Albert Einstein
Center for Fundamental Physics at the University of Bern, and in part by the
Swedish Research Council, project number 2021-06638. A. R.-S. is funded in
part by MIUR contract number 2017L5W2PT.

## References

  * [1] J. Bijnens, N. Hermansson-Truedsson and A. Rodríguez-Sánchez, _Short-distance constraints for the HLbL contribution to the muon anomalous magnetic moment_ , _Phys. Lett._ B798 (2019) 134994 [1908.03331].
  * [2] J. Bijnens, N. Hermansson-Truedsson, L. Laub and A. Rodríguez-Sánchez, _Short-distance HLbL contributions to the muon anomalous magnetic moment beyond perturbation theory_ , _JHEP_ 10 (2020) 203 [2008.13487].
  * [3] J. Bijnens, N. Hermansson-Truedsson, L. Laub and A. Rodríguez-Sánchez, _The two-loop perturbative correction to the $(g-2)_{\mu}$ HLbL at short distances_, _JHEP_ 04 (2021) 240 [2101.09169].
  * [4] T. Aoyama et al., _The anomalous magnetic moment of the muon in the Standard Model_ , _Phys. Rept._ 887 (2020) 1 [2006.04822].
  * [5] Muon g-2 collaboration, _Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL_ , _Phys. Rev. D_ 73 (2006) 072003 [hep-ex/0602035].
  * [6] Muon g-2 collaboration, _Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm_ , _Phys. Rev. Lett._ 126 (2021) 141801 [2104.03281].
  * [7] S. Borsanyi et al., _Leading hadronic contribution to the muon magnetic moment from lattice QCD_ , _Nature_ 593 (2021) 51 [2002.12347].
  * [8] G. Colangelo et al., _Prospects for precise predictions of $a_{\mu}$ in the Standard Model_, 2203.15810.
  * [9] M. Hayakawa and T. Kinoshita, _Pseudoscalar pole terms in the hadronic light by light scattering contribution to muon g - 2_ , _Phys. Rev. D_ 57 (1998) 465 [hep-ph/9708227].
  * [10] M. Hayakawa and T. Kinoshita, _Comment on the sign of the pseudoscalar pole contribution to the muon g-2_ , hep-ph/0112102.
  * [11] J. Bijnens, E. Pallante and J. Prades, _Analysis of the hadronic light by light contributions to the muon g-2_ , _Nucl. Phys. B_ 474 (1996) 379 [hep-ph/9511388].
  * [12] J. Bijnens, E. Pallante and J. Prades, _Comment on the pion pole part of the light by light contribution to the muon g-2_ , _Nucl. Phys. B_ 626 (2002) 410 [hep-ph/0112255].
  * [13] G. Colangelo, M. Hoferichter, M. Procura and P. Stoffer, _Dispersion relation for hadronic light-by-light scattering: theoretical foundations_ , _JHEP_ 09 (2015) 074 [1506.01386].
  * [14] G. Colangelo, M. Hoferichter, M. Procura and P. Stoffer, _Dispersion relation for hadronic light-by-light scattering: two-pion contributions_ , _JHEP_ 04 (2017) 161 [1702.07347].
  * [15] K. Melnikov and A. Vainshtein, _Hadronic light-by-light scattering contribution to the muon anomalous magnetic moment revisited_ , _Phys. Rev._ D70 (2004) 113006 [hep-ph/0312226].
  * [16] P. Masjuan and P. Sánchez-Puertas, _Pseudoscalar-pole contribution to the $(g_{\mu}-2)$: a rational approach_, _Phys. Rev._ D95 (2017) 054026 [1701.05829].
  * [17] M. Hoferichter, B.-L. Hoid, B. Kubis, S. Leupold and S.P. Schneider, _Dispersion relation for hadronic light-by-light scattering: pion pole_ , _JHEP_ 10 (2018) 141 [1808.04823].
  * [18] A. Gérardin, H.B. Meyer and A. Nyffeler, _Lattice calculation of the pion transition form factor with $N_{f}=2+1$ Wilson quarks_, _Phys. Rev._ D100 (2019) 034520 [1903.09471].
  * [19] G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub and P. Stoffer, _Longitudinal short-distance constraints for the hadronic light-by-light contribution to $(g-2)_{\mu}$ with large-$N_{c}$ Regge models_, _JHEP_ 03 (2020) 101 [1910.13432].
  * [20] V. Pauk and M. Vanderhaeghen, _Single meson contributions to the muon‘s anomalous magnetic moment_ , _Eur. Phys. J. C_ 74 (2014) 3008 [1401.0832].
  * [21] I. Danilkin and M. Vanderhaeghen, _Light-by-light scattering sum rules in light of new data_ , _Phys. Rev. D_ 95 (2017) 014019 [1611.04646].
  * [22] F. Jegerlehner, _The Anomalous Magnetic Moment of the Muon_ , vol. 274, Springer, Cham (2017), 10.1007/978-3-319-63577-4.
  * [23] M. Knecht, S. Narison, A. Rabemananjara and D. Rabetiarivony, _Scalar meson contributions to $a_{\mu}$ from hadronic light-by-light scattering_, _Phys. Lett. B_ 787 (2018) 111 [1808.03848].
  * [24] G. Eichmann, C.S. Fischer and R. Williams, _Kaon-box contribution to the anomalous magnetic moment of the muon_ , _Phys. Rev. D_ 101 (2020) 054015 [1910.06795].
  * [25] P. Roig and P. Sanchez-Puertas, _Axial-vector exchange contribution to the hadronic light-by-light piece of the muon anomalous magnetic moment_ , _Phys. Rev. D_ 101 (2020) 074019 [1910.02881].
  * [26] E.-H. Chao, R.J. Hudspith, A. Gérardin, J.R. Green, H.B. Meyer and K. Ottnad, _Hadronic light-by-light contribution to $(g-2)_{\mu}$ from lattice QCD: a complete calculation_, _Eur. Phys. J. C_ 81 (2021) 651 [2104.02632].
  * [27] I.I. Balitsky and A.V. Yung, _Proton and Neutron Magnetic Moments from QCD Sum Rules_ , _Phys. Lett._ 129B (1983) 328.
  * [28] B.L. Ioffe and A.V. Smilga, _Nucleon Magnetic Moments and Magnetic Properties of Vacuum in QCD_ , _Nucl. Phys._ B232 (1984) 109.
  * [29] M.A. Shifman, A.I. Vainshtein and V.I. Zakharov, _QCD and Resonance Physics. Theoretical Foundations_ , _Nucl. Phys._ B147 (1979) 385.
  * [30] K. Melnikov and A. Vainshtein, _On dispersion relations and hadronic light-by-light scattering contribution to the muon anomalous magnetic moment_ , 1911.05874.
  * [31] G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub and P. Stoffer, _Short-distance constraints on hadronic light-by-light scattering in the anomalous magnetic moment of the muon_ , _Phys. Rev. D_ 101 (2020) 051501 [1910.11881].
  * [32] J. Leutgeb and A. Rebhan, _Axial vector transition form factors in holographic QCD and their contribution to the anomalous magnetic moment of the muon_ , _Phys. Rev. D_ 101 (2020) 114015 [1912.01596].
  * [33] J. Lüdtke and M. Procura, _Effects of Longitudinal Short-Distance Constraints on the Hadronic Light-by-Light Contribution to the Muon $g-2$_, _Eur. Phys. J. C_ 80 (2020) 1108 [2006.00007].
  * [34] M. Knecht, _On some short-distance properties of the fourth-rank hadronic vacuum polarization tensor and the anomalous magnetic moment of the muon_ , _JHEP_ 08 (2020) 056 [2005.09929].
  * [35] P. Masjuan, P. Roig and P. Sanchez-Puertas, _A different viewpoint on the Hadronic light-by-light tensor short-distance constraints_ , 2005.11761.
  * [36] L. Cappiello, O. Catà and G. D’Ambrosio, _Scalar resonances in the hadronic light-by-light contribution to the muon (g-2)_ , _Phys. Rev. D_ 105 (2022) 056020 [2110.05962].
  * [37] G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub and P. Stoffer, _Short-distance constraints for the longitudinal component of the hadronic light-by-light amplitude: an update_ , _Eur. Phys. J. C_ 81 (2021) 702 [2106.13222].
  * [38] I. Danilkin, M. Hoferichter and P. Stoffer, _A dispersive estimate of scalar contributions to hadronic light-by-light scattering_ , _Phys. Lett. B_ 820 (2021) 136502 [2105.01666].
  * [39] J. Leutgeb, J. Mager and A. Rebhan, _Holographic QCD and the muon anomalous magnetic moment_ , _Eur. Phys. J. C_ 81 (2021) 1008 [2110.07458].
  * [40] J. Leutgeb and A. Rebhan, _Hadronic light-by-light contribution to the muon g-2 from holographic QCD with massive pions_ , _Phys. Rev. D_ 104 (2021) 094017 [2108.12345].
  * [41] M. Zanke, M. Hoferichter and B. Kubis, _On the transition form factors of the axial-vector resonance f 1(1285) and its decay into e+e-_, _JHEP_ 07 (2021) 106 [2103.09829].
  * [42] J. Aldins, T. Kinoshita, S.J. Brodsky and A.J. Dufner, _Photon - photon scattering contribution to the sixth order magnetic moments of the muon and electron_ , _Phys. Rev._ D1 (1970) 2378.
  * [43] G. Buchalla, A.J. Buras and M.E. Lautenbacher, _Weak decays beyond leading logarithms_ , _Rev. Mod. Phys._ 68 (1996) 1125 [hep-ph/9512380].
  * [44] M. Knecht, S. Peris, M. Perrottet and E. De Rafael, _Electroweak hadronic contributions to the muon (g-2)_ , _JHEP_ 11 (2002) 003 [hep-ph/0205102].
  * [45] A. Vainshtein, _Perturbative and nonperturbative renormalization of anomalous quark triangles_ , _Phys. Lett. B_ 569 (2003) 187 [hep-ph/0212231].
  * [46] A. Czarnecki, W.J. Marciano and A. Vainshtein, _Refinements in electroweak contributions to the muon anomalous magnetic moment_ , _Phys. Rev._ D67 (2003) 073006 [hep-ph/0212229].
  * [47] M. Knecht, S. Peris, M. Perrottet and E. de Rafael, _New nonrenormalization theorems for anomalous three point functions_ , _JHEP_ 03 (2004) 035 [hep-ph/0311100].
  * [48] M. Knecht and A. Nyffeler, _Hadronic light by light corrections to the muon g-2: The Pion pole contribution_ , _Phys. Rev. D_ 65 (2002) 073034 [hep-ph/0111058].
  * [49] M.A. Shifman, _Anomalies and Low-Energy Theorems of Quantum Chromodynamics_ , _Sov. Phys. Usp._ 32 (1989) 289.